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THREE- DIMENSIONAL ELASTIC- PLASTIC FINITE-ELEMENT 
ANALYSIS OF FATIGUE CRACK PROPAGATION 


BY 

R. G. Chermahini 1 and G. L. Goglia 2 
INTRODUCTION 

Fatigue cracks have been a major problem in designing structures sub- 
jected to cyclic loading. Cracks frequently occur in structures such as 
aircraft and spacecraft. The inspection intervals of many aircraft struc- 
tures are based on crack-propagation lives. Therefore, improved prediction 
of propagation lives under flight-load conditions (variable-amplitude load- 
ing) are needed to provide more realistic design criteria for these struc- 
tures. 

The main thrust of this study was to develop a three-dimensional, non- 
linear, elastic-plastic, finite element program capable of extending a 
crack and changing boundary conditions for the model under consideration. 
The finite-element model is composed of 8-noded (linear-strain) isopara- 
metric elements. In the analysis, the material is assumed to be elastic- 
perfectly plastic. The cycle stress-strain curve for the material is shown 
in Fig. 1. Zi enkiewicz' s "initial-stress" method, von Mises's yield criter- 
ion, and Drucker’s normality condition under small-strain assumptions are 
used to account for plasticity. The three-dimensional analysis is capable 
of extending the crack and changing boundary conditions under cyclic 
loading. Initially, the crack is assumed to grow as a straight-through 
crack . 

Research Associate, Department of Mechanical Engineering & Mechanics, Old 
Dominion University, Norfolk, VA 23508. 

2 Eminent Professor, Department of Mechanical Engineering & Mechanics, Old 
Dominion University, Norfolk, VA 23508. 



Using a three-dimensional nonlinear computer program on a cyber-nos 
system was impossible due to its limited storage capacity. To avoid this 
problem, the next alternative was to utilize a VPS- 32 machine with unlimited 
storage capacity. Using the scalar version of the program on the VPS-32 was 
costly due to the plasticity part of the program. Therefore, in order to 
reduce the cost of the computations, the three-dimensional computer program 
was vectorized. 

The finite-element formulation of the program using an 8-noded linear 
isoparametric cubic element is listed in Appendix A. The description of the 
nonlinear program is attached in Appendix B. A list of the program is shown 

N 

in Appendix C. 
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Figure l. CYCLIC STRESS-STRAIN CURVE FOR AN ELASTIC-PERFECTLY 

PLASTIC MATERIAL 
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APPENDIX A 


FINITE-ELEMENT FORMULATION 


il 




Figure 2. Arbitrary hexahedron. 

(a) Elastic Analysis 

The basic concept of the finite-element method is that any continuous 
quantity can be approximated by a discrete model composed of a set of 
piecewise continuous functions defined over a finite number of subdomains 
( 1 ). 

An isoparametric 8-noded cubic element (Fig. 2.) was utilized in the 
formulation of the elastic-plastic structure into the nonlinear computer 
program. The following section describes cubic element. 

Displacement Functions . The displacement function (2) for any point in the 
cubic element is defined as: 
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u(x,y,z) = a x + a 2 x + a 3 y + a 4 z + a 5 xy + a 6 yz + a 7 xz + 'a 8 xyz 

v(x,y,z) = b 2 + b 2 x + b 3 y + b 4 z + b 5 xy + xy + b 6 yz + b 7 xz + b 8 xyz 

w(x,y,z) = c 2 + c 2 x c 3 y + c 4 z + c 5 xy + c 6 yz + c 7 xz + c 8 xyz A(l) 

where u, v and w are displacement in the x, y and z directions, 
respectively. The constant coefficients are determined by imposing the 
nodal coordinates of each cubic element into equations A (1). The above 
displacement function can be applied to the cubic element as long as the 
sides of the cubic element are defined by planes parallel to the coordinate 
planes. However, for the elements whose sides are skewed, the above dis- 
placement function no longer is applicable. Therefore, in order to avoid 
this restriction, an 8-noded linear isoparametric cubic element is employed 
(Fig. 2.). 

The original cube can be mapped on to a cube of 2x2x3 unit (2) in the 
?, n, K space by the transformation 

x = aj + a 2 s + a 3 n + a 4 s + a 5 cn + a6hC + azCS + as? 1 !? 

y = bi + b 2? + b 3T) + b 4? + b 5?T1 + b 6ri £ + b 7Z e, + b 8?T1 s A(2) 

z = ci + c 2 ? + c 3 n + c 4 £ + c 5 5n + c 6 n? + c 7 z,l + cs^S 

The values of the coefficients in equation A(2) depend on the nodal coordin- 

ates of each cubic element and are different for different elements. The 
transformation is defined by polynomials in c, n and 5 which is contin- 
uous within the element, the continuum confined within an element in x, y 
and z coordinates is mapped on to a continuum within the 2x2x2 cube in 5, n 
and % coordinates. It remains to be shown that the transformation is 
continuous across two adjoined elements, that a common surface between 
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two adjoined elements in the x, y, z space will transform into a common 
surface of two adjoined cubes in n, 5 space. 

If we assign the following values of the parameters n, ? to 
the faces of distorted elements shown in (Fig. 3.) one yields: 



Face 

Coordinate value 

pokl 

S = 1 

mnji 

c = 1 

impl 

n = 1 

jnok 

n = -1 

mnop 

5 = 1 

i jkl 

e = i 


Fig. 3. Linear isoparametric cubic element. 

Therefore, the nodal points i , j , k , 1 and m, n, o, p will have the 
following coordinates in the ?, n, 5: 
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nodal point 



Now the displacements (u, v, w) in 
as: 


coordinates 



x, y, z directions can be written 


11 = 0! + a 2 c + a 3 n + 45+ a s 5n + a 6 nS + a 7 £5 + a 8 ?nC 
v = $! + g 2 ? + s 3 n + 64? + fs?o + MS + M? + 2 8 ?n? A(3) 

w = yi + Y2? + Y 3 n + Y 4? + Ys? + + Y7?? + Ye?^? 

which are continuous (2) within the elements as well as across the surfaces 
common to any two adjoined elements. Consider the term u in equations 
A(3), denote by {a} and {U} the vectors for the a's and u^ nodal 
displacements of all the nodal points of the element. Inserting the values 
of u • , and 5^ for the various nodal points, we obtain eight 

equations corresponding to the first equation of A(3) which can be written 
as 


{U} - [Ai] {a} 


A ( 4 ) 
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Let's define [o^J = 

[Aj], and 

thus 

have 

{*} " 

L ai J {u}. Now the 

displacement functions 

for the distorted i 

element 

can be 

written as: 


u- LSJ 

C«i] 

M 




v= LSJ 

W 



A(5) 


w= lsj 

[«i] 

(wl 




where [S] is defined as: 

LSJ = La C n 5 ?n ?5 ?n?J A(6) 

The shape functions for the isoparametric 8-noded element can be determined 
(1) from the product of [S] and [aj] matrices. 

N i =1 (l+ ?i ) (l+n i ) (1+Ci) A(7) 

8 

where n^, C i = ± 1 and i = 1, 2, 8. 

The x, y, and z coordinates at any point in the element, can be 
expressed in terms of shape functions : 

8 

x = y n. x. 

'll 

1=1 

8 

y = l N. y. A(8) 

i=l 1 1 

■ * ■ Li 2 < 
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’ , 

X 


[N] 

Co] 

CO] 



y 

' = 

[0] 

[N] 

CO] 


fy„) 

z 


[0] 

CO] 

CN] 


$ 


where {x n } T » [ Xl x 2 ... x 8 j, {y n } T - [y x y 2 ... y 8 j 


and 




Czi z 2 


Zfl]' 


Element Strain ; The elastic strain at any point within the element is given 
by [3] 


{e} 


f e 1 


r 3U > 

! 


3x 

) 


3v 

E y 


3y 

£ 


3w 

z 


— 


► < 

32 , 

Y 


3u 3v 

xy 

= 




3y 3x . 

Y 


3v , 3w 

yz 


+ 



3z 3y 

^ Y zx •> 


3w , 3u 
v. — + — J 


v — f — j 

3x 3z 


[ Bj ] [B 2 ].-..[BgJ 


{u} = [Bj {u} 


A(9) 


where the matrix [B] is defined as: 
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A(10) 


The transformation relationship between local and global coordinates is 


given by: 



where [J] is the Jacobian matrix and it is defined as: 



where { x n } T= C x i x 2 ••• x 8^* 


A(ll) 


A (1 2) 
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Element Stress. For linear-elastic and isotropic materials, the element 
stresses are calculated using Hook's law 


{a} - [D] { e } + fd} 


A( 1 3) 


The strain vector is { e} = [B] {u}; therefore, the stresses are 


{a} = [D] [B] {u} + {d} 


A( 14) 


where {a } is initial stress which may exist in the element. The material 
property matrix [D], is defined as: 


[ 0 ] = 


(I4v)(l-2v) 


1 -v v v 0 0 0 


1 — V V 0 0 0 


Symmetrical 


1 -v ooo 

1 -2v o o 


l-2v 


l-2v 


A(15) 


where E is Young's modulus and v is poisson's ratio for the material. 
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Element Equation s. The potential energy Tr p (u, v, w) which is composed of 
strain energy Up (u, v, w) and v p (u 4 v, w) the work done by the applied 
loads during displacement changes is given by [3]. 

* * - ///„ M T M {c}(iv - Ilf [F*] (5} dv + If [T*] (7) ds A(16) 

p 2 v v b i 

★ ★ * * ★ ★★★ 
where [F ] = [x y z ], [T ] = [T x T y T z ]. 

and [s] = [u,v,w]. 

The equilibrium equations for the element are obtained by taking partial 
derivatives of * with respect to u x , v x , w x , etc., and equating to 
zero. 


EL = 

3{U} 


0 


> 


which leads to the 24 element equilibrium equations as 


[K] { u} ~ {q} = { Qi} + {q 2 l 

24x24 24x1 24x1 24x1 24x1 


where [K] is the element stiffness matrix, 

M ■ III, [B] T [0] [B] dv + [K s ] 

and { Q} is the element nodal load vector. 


A( 1 7) 


A(18) 


A(19) 
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tQ} - IQi} + {q 2 } - III, [n] t {F*} <t* + // Si [N] T (T*} ds 


A( 20 ) 


The diagonal matrix [K $ ] in Eq . A(19) is the eleastic stiffness of the 
springs, which are connected to the boundary nodes. 

(b) Elastic-plastic analysis 

Finite-element techniques applied to linear elastic materials have been 
solved successfully. However, for an elastic-plastic material, the coeffi- 
cient in the stiffness matrix varies as a function of material loading. Two 
computational methods have been used successfully in the solution of 
elastic-plastic problems. In the first, the change at each step of load 
increase in plastic strain is calculated and treated as an initial strain 
for which the elastic stress distribution is adjusted (1). This method 
fails if ideal plastic is postulated or if the degree of hardening is small. 
In the second method, the "incremental stress method," the stress-strain 
relationship for every load increment is adjusted to account for plastic 
deformations. The work of Pope (4), Swedlow (5), Marcal and King (6), Reyes 
and Deere (7) and Popov and others (8) falls into this category. 

The "incremental elasticity" method has one serious disadvantage. At 
each step of the computation the stiffness matrix of the structure is up- 
dated and iterative schemes of solution are necessary to avoid excessive 
computational costs. To minimize computational costs, the "initial stress" 
approach is used (1). In the incremental stress method, the basic elasti- 
city matrix remains unchanged. This technique converges more rapidly than 
the initial strain method. 

Yield Criterion . In any elastic-plastic analysis, it is necessary to intro- 
duce a yield criterion to determine the state of stress at which yielding 
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r 

occurs. The von Mi ses yield criterion or maximum distortion energy theory 

Is 

of failure, which finds considerable experimental support in ductile mater- 
ials, is used to determine whether the material at any point in the 
structure has yielded. This criterion assumes that yielding begins when the 
distortion energy equals the distortion energy at yield in simple tension 
(1). The von Mi ses yield criterion for a three dimensional state of stress 
is given by 


F = F ( CT ) = [I (a x -cr y ) 2 + j (<r y -a z ) 2 + \ ^Z^X^ 

l 

+ 3T 2 + 3T 2 + 3T 2 1 2 - a A(21) 

xy xz xz J 

where ~o = o’ (K) is the uniaxial stress at yield. If F (a) < 0, the 
material is in elastic range. If F (a) > 0, the material has experienced 
plastic deformation and one of the flow theories of plasticity must be used 
for determining the components of plastic strains and stresses due to the 
applied load. 

During an infinitesimal increment of stress, changes of strain are 
assumed to be divisible- into elastic and plastic parts (1). Thus, the 
strain increment can be written as: 

{de} = {de g } + {de p } A(22) 


where the elastic strain increments are related to the stress increments by 
the symmetric material matrix D. The plastic strain increments are related 
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to the yield criterion through Drucker's normality principle 


{■V = x f-f 

! Scr 


A(23) 


Therefore; Eq. A(22) can be rewritten as: 


Ids} = [D I - ' ftto} + X { 11 } A(24) 

3o 


At the point of incipient plasticity, the stresses are on the yield surface 
and the yield function is given by: 


F(0,k) = 0 


A(25) 


where K is a hardening parameter. 

Differentiating A(25) results in: 

d = II 5a! + d a2 + . . . + 11 dK = 0 A(26) 

9cT l 302 3 k 

or {—} Tdo - AA = 0 . A(27) 

30 

Solving for A gives 


9F A If 1 

= dk - 


3k A 


A(28) 


15 



Equations A(24) and A(27) can be written in matrix form as 


de 1 


d" 1 

9F 



1 

. 0 , 

\ 

= 

(-) T 
9a 

9a 

-A 

1. ) 

A(29) 


The constant x can be eliminated from Eq. A(23). The final expression 

which relates the stress changes in terms of imposed strain changes can be 

* 

written as: da = de 

ep 

A(30) 

or 


where 


and 


* 


D-0(— ! D 


A + {— } D 1 


J-l 


9 a 


9 a 


9 a 


9 a 


$ “ [F X F y F z F x, F y Z F xzJ 


A + 


3a i 

3a 2 

3a 3 

= , 

F y = — ’ 

F z = 

2 a 

2 a 

2 a 


3T 

3T 

3T 

xy 

F = yZ 

F = ZX 

» 

yz 

’ zx 

"a 

a 

o’ 


in which the dashes stand for deviatoric stresses i.e. 


A (31) 


A(32) 
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u 

★ 

The elastic-plastic matrix 0 g p replaces the elastic matrix D in 
incremental elastic-plastic analysis. The plastic load vector for the 
elements which deform plastically is given by: 

{dq} = fff CB3 T {da} dv m A(33) 

where {da} is defined as: 

{do}= {do g } - { dcr } + ( [De] - [Dep]) { de} A(34) 
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APPENDIX B 


Description of the Finite-Element Computer Program 
The computer program presented here was based on the three-dimensional 
8-noded linear isoparametric cubic element. The optimum goal of this study 
was to develop a three-dimensional nonlinear computer program capable of 
extending a crack and changing the boundary conditions for the model under 
consideration. This program in its present form is not a general analysis 
program for nonlinear cracked structures. The restrictions are listed as 
follows: (1) the crack must lie on the x-axis and propagate in the positive 

x-direction, (2) the configuration and loading must be symmetric about tne 

x-axis. 

The input to the program is illustrated by using one eighth of a 
center-crack panel shown in Fig. 4. 

1. CRACK, WIDTH, THICK, HEIGHT, DAX:, SCALE (6E10.4) 

The format for each input is shown in parenthesis. Crack specifies the 
crack length in the y=0 plane. Width, thick, height represent width, 
thickness and height of the structure. , DAX is defined as the smallest 
element size in the region and is used for the crack-extension in the 

program. Scale, scales the width, thickness and height of the specimen 

to the desired dimension. 

2. LPRIT, LMAX, KMAX, NLAYER, NEP (1615) 

LPRIT = 0 indicates that no intermediate output is printed. LPRIT = 1 
results in intermediate output. LMAX is the number of nodes in Z=0 
plane. KMAX is the number not elements in Z=0 plane. NLAYER indicates 
the number of layers in the structure. NEP specifies elastic or 
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plastic anaslysis if NEP=0, elastic analysis is performed. If NEP 
> 0, the plastic analysis is performed. 

3. K, XR ( I ) , YR ( I ) , ZR ( I ) , (15, 4X , 3 E15.7) 

K refers to the node number, and XR(I), YR(I), ZR( I ) are the 
coordinates of node K in x,y and z direction, respectively. 

4. IN, (M0DE(J, IN), J = 1,8) (1615) 

IN describes the element number, and node gives the nodal 
connective of each cubic element in the structure. 

5. NSYMPL (1615) 

NSYMPL specifies the nunber of synmetric planes 

6. (ISYMPLY(I), I = 1, NYSMPL) (1615) 

ISYMPL describes the corresponding nunbers designated for each 
plane in the structure. 

7. NFIX, NLOAD, SNPD (1615) 

NFIX, NLOAD NSPD describe the number of fixed loaded, and specified 
displacements for nodes, respectively. 

8. NODF , MU, MV, MW (1615) 

NODF describes the number of fixed nodes, and MU, MV and MW 
represents the u, v and w displacements fixed for each node. 

9. Nodlod (IL), P ¥ , P„ P 7 (1615) 

Nodlod specifies the number of loaded nodes, and p , p and p 

x y l 

represent the components of loading in x, y, and z direction, 
respectively. 

10. NODS, K, DISP(N) (1615) 

NODS is the node number, K is the code for u, v and w 
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displacements, and Disp is the specified displacement for the 
corresponding node. 

11. NTYP, NLM, SCRIT, RP, ACURCY (215, 4E10.4) 

NTYP stands for the crack growth criterion. NLM is the number of 
increments to release the crack tip force. SCRIT is used for the 
CTOD criterion. RP is the relaxation parameter and ACURCY is used 
for the crack opening displacement accuracy. 

12. P, WORD (E103, IX, \) 

P designates the maximum applied stress for each cycle. The word 
specifies stationary or growing crack for each cycle. If word is 
set equal to grow, the crack will extend one element size. If word 
is equal to halt, the crack will be stationary for that cycle. 
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APPENDIX C 


FORTRAN LISTING 



on o oononnoonnnoonnnnooooonnnoo oonooononnn 


ORIGINAL PAGE SS 

OF POOR QUALITY 


PROGRAM CRACK1 ( INPUT , OUTPUT , TAPE7-DI , TAPE5-INPUT , 

1TAPE6-0UTPUT) 

. COMMON/ MAIN/ AA( 2000000) ,BB<9600,1) ,D(6,-6) ,DINV(6,6) , 

1DISP(80) ,EPS(1240C) ,EFEST(I550) ,F0RCE(10) ,LINE(100) , 

2L0CAT( 10) .LBFOR(.Vj) ,MB(9600) ,MSUM(9600) ,MPTAB(9600) , 

• 3HPLAS/ 12400) ,M0DE(8 ,1550) ,MPLC(1550) ,NODXO(700) , 

4NODYO(700) ,NODZO(700) ,N0DXC(700) ,N0DYC(700) ,N0DZC(700) , 

5N0DFIX(80) ,N0DL0D(80) ,NDISP(80) ,R(9600) , 

6SIGBAR( 12400) ,SK(24,24) ,T1<9600) ,T2(9600) ,T4(2000) , 

7T3(9600) ,U(3200) ,U0LD(3200) ,V(3200) ,VOLD(3200) ,V2(3200) , 

8W(3200) ,WOLD(32O0) ,X(74400) ,XR(3200) ,Y(74400) , 

9YR(3200) ,Z(9600) ,ZR(3200) 

COMMON/ CNST/ EPSI , SK2 , LMAX, KMAX, DAX, 

1PYLE, SCRIT , YOUNG , POIS , CRACK, PT , WIDTH , PMAX , HP , 

2SBAR, LPRIT , NGAUS , NLAYER , NNODE , 

3IN0DX0 , INODYO, INODZO , INODXC , INODYC , INODZC ,LNSTIF,MXNOD , 

4MXNEL, MXGAUS , ICUT , LTOTB , ITNODX ,KLU , NTYP , NLM , 

5NDOF , KNEW , NE P , ERIT , NELM , AM , ROM 

COMMON/ MLTNMAT/ YSTRS (20) ,YSTRH(20) ,PU'10DR(20) , NSEGMT 
C0MM0N/D382/NNPE,NDF,NQD,NSTR,HQD2,HNPE2,NQD2NPE,NQD2SR,MXQ2S 
COMMON/ VECT/ STRV( 8 , 6) , STRSV( 8 , 6) , BMT( 64 , 3 ) , WDUM( 8 ) , XE ( 8 , 3 ) , 

C NCUBE(8) ,DIS(S,3) 

DIMENSION IIMAX(3200) ,NSAM£(3200,20) ,MS(8) 

DIMENSION JNEW(3200) ,TITLE(20) ,ISYMPL(6) ,NBEGIN(8) ,NEND(8) 

DIMENSION STR(6 ) 

*********************************** ****** ********** 

* XR(I,J) ,YR(I,J) COORDINATES OF RECTANGULAR ELEMENTS* 

* WHICH ARE LOCATED IN THE Z-0 PLANE. 

* XR(I)',YR(I),ZR(I) COORDINATES OF NODES IN THE STR* 

*UCTURE. * 

*NODXO(I) NODE NUMBERS FOR PLANE X-0 * 

*NODYO(I) Y-0 * 

*NODZO(I) Z-0 * 

*NODXC(I) X-XCOR * 

*NODYC( I) Y-YCOR * 

*NODZC(I) Z-ZCOR * 

*U(I),V(I),W(I) DISPACEMENT COMPONENTS "OR EACH NODE 1/ THE SPECIMEN 
N0DL0D(80) MAX OF 80 NODES LOADED 

* 

EPSI IS ACCURACY CHECK VALUE 

SK2 STIFFNESS OF SPRINGS CONNECTED TO BOUNDARY NODES 

LMAX NO OF NODES IN Z-0 PLANE 

KMAX NO OF ELEMENTS IN Z-0 PLANE 

DAX SMALLEST ELEMENT SIZE IN THE STRUCTURE 

PUD LOAD AT INITIAL YIELD 

SCRIT USED FOR CTOD CRITERION 

YOUNG YOUNGS MODULUS OF THE MATERIAL 

POIS POISSON RATIO OF THE MATERIAL 

CRACK CRACK LENGTH 

PT VARIABLE USED FOR LOADING 

WIDTH WIDTH OF THE SPECIMEN 

SIGYS YIELD STRESS OF THE MATERIAL 

LPRIT LPRIT GREATER THAN 0 NO INTERNAL OUTPUT .LPRIT-0 
INTERNAL OUTPUT(USED FOR SMALL PROBLEMS) 

NGAUS NO GAUSS POINTS IN EACH DIRECTION 
NLAYER NO OF LAYERS PUT IN THE STRUCTURE 
NNODE TOTAL NO OF NODES IN THE STRUCTURE 
INODXO TOTAL NO OF NODES IN X-0 PLANE 

INODYO TOTAL NO OF NODES IN Y-0 PLANE 

■ INODZO TOTAL NO OF NODES IN Z-0 PLANE ' ! 

INODXC TOTAL NO OF NODES IN X-C PLANE 

INODYC TOTAL NO OF NODES IN Y-C PLANE - ! 

' INODZC TOTAL HO OF NODES IN Z-C PLANE 

. LNSTIF MAXIMUM DIMENSION FOR AA MATRIX 

MXNOD MAXIMUM NODES PUT INTO THE PROGRAM 
MXHEL MAXIMUM ELEMENTS PUT INTO THE PROGRAM 

j 
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OR5G5MAL PAGE S3 

OF POOR QUALITY 


C MXiJOD AND HX11EL ARE FOR DIMENSIONAL PURPOSES 
C MXGAUS MAXIMUM NO OF ELEMENTS MULTIPLY THE NO OF GAUSS 
C POINTS IN EACH DIRECTION(X, Y,Z IF NGAUS-2.THE NO IS 2*2*2) 

C ICUI VARIABLE USED IN BREAK SUBPROGRAM FOR RELEASING FORCES 

C LTOTB TOTAL NO NODES IN THE THICNESS ALONG THE CRACK TIP 

C . ITNODX TOTAL NO OF NODES ALONG THE CRACK. LINE 
C KLU VARIABLE USED FOR CRACK EXTENSION 

C NTYP VARIABLE USED FOR TYPE OF CRACK EXTENSION 

C NLM NO OF INCREMENTS TO RELEASE THE NODAL FORCES 
C NLOAD NO OF LOADED NODES IN THE STRUCTURE 

C NSPD NO OF SPECIFIED DISPLACEMENTS 

C MAXIT MAX NO OF ITERATION USED FOR CONVERGENCE PURPOSES 
C NDOF TOTAL NO OF DEGRES OF FREEDOM IN THE MODEL 

C NEP IF NEP -0 ELASTIC ANALYSIS, IF NEP GREATER 0 PLASTIC ANAL 

C ERIT ACCURACY CHECK VALUE FOR CONVERGENCE USED IN SUB PLAS 

C NELM TOTAL NO OF ELEMENTS IN THE SYSTEM 

C AM, ROM LINEAR OR NONLINEAR STRAIN HARDENING COEFFICIENTS 
C IF AM-0 MATERIAL IS ELASTIC-PERFECTLY . 

C KNEW VARIABLE USED IN CONTACT SUBPROGRAM TO CHECK WHETHER 
C THE NODE CLOSED OR OPENED. 

C 

DATA NNPE , NDF , NQD , NSTR/8 ,3,2,6/ 

C *** OPEN MAP AND ZERO THE AA VECTOR OF LENGTH LENTOT 
LENTOT-2000000+9600*9+1550*130+700*(6)+80*4+100 
1+10*3+72+2000+3200*10+576 
J-LENTOT/65536 

JJ-LENTOT-(LENTOT/65536) *65536 

IF(JJ.NE.O) J-J+l NEW 

L0PN-J*128 NEW 

CALL OPEN(LOPN) 

C *** ZEROING THE VECTORS 
J-LENTOT/65536 
DO 223 I-l.J 
I l-(I-l) *65535+1 
AA( II; 65 53 5) -0.0 
223 CONTINUE 

J-J*65536+l 
J J-LENTOT- J+l 
AA(J ; JJ)« 0.0 
C 

cc *** 

c 

LNSTIF-2000000 
MXNEL-1550 
MXN0D-3200 
NGAUS-2 
NQD2-NGAUS**3 
NNPE2-( NNPE* ( NNPE+1 ) ) / 2 
NQD2NPE-NQD2*NNPE 
NQD2SR-NQD2*NSTR 
MXQ2S-MXNEL*NQD2SR 
MXGAUS-NQ2*MXNEL 
LL”3*MXN0D 
MPTAB ( 1 ; LL) -0 
Z(l;LL)-0.0 
R(1;LL)»0.0 
BB( 1 , l;LL)-0 . 0 
ZR( 1 ;MXNOD)-0.0 
CALL Q3CL0CKS( CPU, WALL) 

42 • FORMAT ( 15, 3F1 0.3) 

C 

C **’* READ GEOMETRIC DATA 
C 

READ (5,222) (TITLE(I) ,1-1,20) 

222 FORMAT ( 2 0A4) 

WRITE(6 , 15) (TITLE(I), 1-1,20) 


NEW 
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15 F0RMAT(1H1///5X,20A4) 

READ(5 ,16) CRACK , WIDTH, THICK , HEIGHT , DAX , SCALE 
. WRITE(6 , 17) CRACK , WIDTH , THICK .HEIGHT , DAX .SCALE 

16 . F0RMAT(6E10.4) 

17 FORMAT(5X,'CRACK-',F10.4,2X,'WIDTH-',F10.4,2X,'THICK-',F10.4, 

• C/ / 5X , 'HEIGHT-' , FIO . 4 ,2X,'DAX-' , F10 . 6 , 2X, ' SCALE-' , F10 . 5) 

39 F0RMAT(16I5) 

XCOR-UIDTH 

YCOR-HEIGHT 

ZCOR-THICK 

EPSI-l.E-10 

READ(5 ,39) LPR IT , LMAX , KMAX , NLAYER , MEP 
WRITE(6,28) LPRIT , LMAX, KHAX, NLAYER, NEP 
2 ? FORMAT(5X,'LPR-',I2,2X,'LMAX-',I5,2X,'KMAX-',I5,2X, 

C 'NLAYER-' ,I2,2X,'NEP-' ,12) 

NNODE-( NLAYER+1 ) *LMAX 

NDOF-NNODE*3 

NELK-KMAX*NLAYER 

C CONSTANTS IN POLYNOMIAL AND D-MATRIX 

CALL ACAL 

READ (5,39) NMAT.NSEGHT 
DO 3 I-l.NMAT 

READ(5,16) YOUNG, POIS.SIGYS, AM, ROM 
WRITE (6 , 4 ) YOUNG , POIS.SIGYS, AM, ROM 
READ(5 ,39) (NBEGIN(IG) ,NEND(IG) ,IG-1 ,8) 

WRITE(6,39)(NBEGIN(IG) ,NEND(IG) ,IG-1,8) 

DO 5 IG-1,8 

IF(NBEGINCIG).EQ.O) GOTO 3 

11- NBEGIN(IG)*8-7 

12- NEND(IG) *8-11+1 

5 SIGBAR(I1;I2)-SIGYS 

3 CONTINUE 

4 F0RMAT(//10X, 'MODULUS, NUE, YIELD STRESS, AM, <* SOM: ' .5E12.4) 
CALL DCON(YOUNG,POIS,D,DINV) 

C 

C *** READ COORDINATES AND CONNECTIVITY 

C 

DO 30 I-l.NNODE 
JNEW(I)»I 

30 READ(7 ,20) K,XR(I) ,YR(I) ,ZR(I) 

20 FORMAT(I5,4X,3E15.7) 

WRITE(6 ,333) 

WRITE(6,861)(J,XR(J) ,YR(J) ,ZR(J) ,J«l,NNODE) 

861 FORMAT(2(3X,I5,3(E13.6 ,1X) )) 

333 F0RMAT(1H1//10X, 'NODAL COORDINATES .NODE# , X,Y,AND, Z'//) 

DO 31 IE-l.NELM, 

31 READ(7 ,39) IN,(MODE(J,IN) ,J-1,8) 

WRITE(6,334) 

334 F0RMAT(1H1//5X, 'NODAL CONNECTIVITY IE, I,J,K,L, I1,J1,K1,L1'//) 
WRITE(6,864) ( IE , (M0DE( J ,1E) , J-l ,8) .IE-l.NELM) 

864 F0RMAT(2(5X,9I5) ) 

C 

C *** 

c 

IZIP1-5 

CALL Q3CLOCKS( CPU, WALL) 

WRITE(6,9999) IZIP 1 , CPU , WALL 

9999 F0RMAT(5X, 'STEP#', 13, 2X, 'TIME IN SECS: CPU-' ,F10,4,2X, 

C'WALL-' ,F12.3) 

WRITE(6 , 1607) NELM 

1607. F0R!1AT(5X, 'TOTAL NO OF HEXAHEDRAN-' ,16) 

C *** IDENTIFY NODES ON CONSTANTS PLANES 
C ' IDENTIFY X-0 PLANE .STORE NODXO ARRAY 

C IDENTIFY Y-0 PLANE .STORE NODYO ARRAY 

C IDENTIFY Z-0 PLANE .STORE NODZO ARRAY 

C IDENTIFY X-C PLANE .STORE tiQDXC ARRAY 
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C IDENTIFY Y-C PLANE j STORE HODYC ARRAY 

C IDENTIFY Z-C PLANE .STORE NODZC ARRAY 

C 

ItjODXO-O 
INODYO-O 
. 'INODZO-O 
INODXC-O 
INODYOO 
INODZC-O 

DO 1300 I-l.NNODE 
IF(ABS(XR(I)) .LE.EPSI) GO TO 1301 
DX1-ABS( XR( I )“XCOR) 

IF(DXl.GT.EPSI) GO TO 1302 
INODXC-INODXC+1 
NODXC( INODXC) “I 
GO TO 1302 

1301 INODXO-INODXO+l 
NODXO(INODXO)-I 

1302 IF(ABS(YR(1) ) .LE.EPSI) GO TO 1303 
DY1»ABS(YR(I)-YC0R) 

IF(DYl.GT.EPSI) GO TO 1304 

INODYOINODYC+1 
N0DYC(IN0DYC)-I 
GO TO 1304 

1303 INODYO- INODYO+l 
NODYO(INODYO)-I 

1304 IF(ABS(ZR(I)). LE.EPSI) GO TO 1305 
DZ1-ABS( ZR( I )— ZCOR) 

IF(DZl.GT.EPSI) GO TO 1300 

INODZC-INODZC+1 


1305 


1300 


1002 


1122 


1001 


NODZC(INODZC)-I 
GO TO 1300 

INODZO-INODZO+1 
NODZO ( INODZO) “I 
CONTINUE 
WRITE(6 , 1002) 

FORMAT( 5X , ' INODXO , 5X, INODYO , 5X, INODZO , 5X, INODXC , 5X, INODYC 
1 , 5X.IN0DZC') 

WRITE ( 6,1122) INODXO , INODYO , INODZO , INODXC , INODYC , INODZ C 
F0RMAT(8X,6I6) 

(NODXO(I).I-l, INODXO) 

(NODYO(I) ,1-1, INODYO) 
(NODZO(I).I-l, INODZO) 

(NODXC(I) ,1-1, INODXC) 
(NODYC(I).I-l, INODYC) 

(NODZCCI) ,1-1 .INODZC) 


IF(LPRIT.EQ.O) 

IF(LPRIT.EQ.O) 

IF(LPRIT.EQ.O) 

IF(LPRIT.EQ.O) 

IF(LPRIT.EQ.O) 

IF(LPRIT.EQ.O) 

CONTINUE 
IZIP1-10 

CALL Q3CL0CKS( CPU, WALL) 
WRITE(6,9999) IZIP1, CPU, WALL 


WRITE(6 ,39) 
WRITE(6,39) 
WRITE(6,39) 
WRITE(6,39) 
WRITE(6 ,39) 
WRITE(6 ,39) 


C 

C **** 

c 


CALL NS AMC ( NODE , NSAME , UMAX , !1B , NNODE , NELM , 8 , MXNOD , MXNEL , NDOF ) 

• MSUH(l)-0 
MSUM(2)-1 
DO 352 1-3, NDOF 
LN-I-1 

352 MSUM( I )-MSUM( LN)+MB(LN) 

LDOF-tlSUM(NDOF)+MB(NDOF) 

WRITE(6 ,504) LDOF 

504 F0RMAT(//10X, 'STORAGE REQUIREMENT FOR STIFFNESS MATRIX IS-'IIO) 

‘ IZIP1-14 


CALL Q3CL0CKS( CPU, WALL) 
WRITE(6,9999) IZIP1, CPU, WALL 
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c 

C ASSEMBLE THE STIFFNESS MATRIX !C 

C 

DO 943 I-l.NELM 
. 'NCUBEU;8)-MODE(l,I;8) 
MS(l;d)-NCU3E(l;8) 

CALL CORDIN(NCUBE,MXNOD,XR,YR,ZR,XE) 
CALL SMALLK(SK,XE,D,IERR) 

DO 943 J»L,8 
DO 943 L-i,8 

IF(MS(L) .LT.MS(J)) GO TO 943 
IU-3*MS(J)-2 
IV-IUfl 
IW-IV+1 
JU-3*MS(L)-2 
JV-JIH-l 
JH-JV+1 

N1-MSUM(JU)-JU+MB(JU)+IU 

N2-N1+1 

N3-N2+1 

N4-MSUM(JV)-JV+MB(JV)+IU 

N5-N4+1 

N6-N5+1 

N7-MSUM(JW)-JW+MB(JW)+IU 

N8-N7+1 

N9-N8+1 

MCl-3*J-2 

MC2-MC1+1 

MC3-MC2+1 

MRl-3*L-2 

MR2-MR1+1 

MR3-MR2+1 

AA(H1)-AA(N1)+SK(MR1,MC1) 

AA( N4) - AA( N4 )+SK( MR2 , MCI ) 

AA( NS ) -AA( N5 )+SK( 11R2 , MC2 ) 

AA( N7) -AA( N7 )+SK( MR3 ,MC1 ) 

AA( N8 ) - AA( N8 )+SK( MR3 , MC2 ) 
AA(N9)-AA(N9)+SK(HR3,MC3) 

952 IF(J.EQ.L) GO TO 943 

AA( N2 ) - AA( N2)+.TK( MR1.MC2) 

AA( N3 ) - AA( N3 )+SK( HR1 , MC3 ) 

AA( N6 ) -AA( N6 )+SK( MR2 , MC3 ) 

943 CONTINUE 

IZIP1-18 

CALL Q3CL0CKS(CPU HALL) 

HRITE(6,9999) IZIP1, CPU, HALL 
C 

C *** IMPOSE SYMMETRIC BOUNDARY CONDITIONS 
C 


315 


316 


READ (5,39) NSYMPL 
WRITE(b,315) NSYMPL 

FORMAT(/5X.,' If OF SYMMETRIC BOUNDARY CONDITIONS -',13) 
IF(NSYMPL.EQ.O) GOTO 314 
READ(5,39) (ISYMPL(I) ,1-1, NSYMPL) 

WRITE (6, 31 6) (ISYMPL(I) ,1-1, NSYMPL) 

FORMAT(10X,' SYMMETRIC PLANE NUMBERS ARE :',6I3) 

DO 317 IS-1, NSYMPL 
ISY-ISYMPL(IS) 

IF(ISY.EQ.l) CALL SY11PLN( AA,MSUM,MB ,MPTAB ,NODXO, INODXO , SK2 , 1 , 
CNDOF.LNSTIF) 

IF(ISY.EQ.2) CALL SYMPLN(AA,MSUM,MB,MPTAB,NODYO,INODYO,SK2,2, 
' CNDOF.LNSTIF) 

IF(ISY.EQ.3) CALL SYMPLN(AA,MSUM,MB,MPTAB,NODZO,INODZO,SK2,3, 
CNDOF.LNSTIF) 

IF(ISY.EQ.4) CALL SYMPLN( AA,MSUM,MB ,MPTAB ,N0DXC , INODXC , SK2 , 1 , 
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CNDOF.LNSTIF) 

IF(ISY.EQ.5) CALL SYMPLN( AA, MSUM , MB , MPT AB , NODYC , INODYC , SK2 , 2 , 
, CNDOF.LNSTIF) 

IF(,ISY.EQ.6) CALL SYMPLN(AA,MSUM,MB,1IPTAB,N0DZC,IN0DZC,SK2,3, 
CNDOF, LNSTIF) 

317 • CONTINUE 
314 CONTINUE 
C 

C *** SYMMETRIC BOUNDARY CONDITIONS ON THE CRACK PLANE 
C 

DO 318 I-l.INODYO 
L-NODYO(I) 

SAP-XR(L) 

IF(SAP.LT. CRACK) GO TO 318 
NV-3*L-1 

NVNV-MSUII(NV)+MB(NV) 

MPTAB(NV)-MV 

AA( NVNV ) - AA( NVNV )+SK2 

318 CONTINUE 
C 

IZIP1-27 

CALL Q3CLOCKS( CPU .WALL) 

WRITE(6,9999) IZ1P1, CPU, WALL 
C ***** READ BOUNDRAY CONDITIONS AND LOADING 
C 

C *** FIXED NODES AND LOADING 
C 

READ(5 ,39) NFIX , NLOAD , NSPD 
WRITE (6, 40) NFIX, NLOAD, NSPD 

40 F0RMAT(//5X,'// OF NODES: FIXED-' ,13 ,2X, 'LOADED-' ,13 ,2X, 

C 'SP. DISP-' ,13//) 

IF(NFIX.EQ.O) GOTO 417 
DO 416 IFIX-1 ,NFIX 
READ(5 ,39)HODF,ttU,MV,MW 
WRITE (6,39) tlODF , MU , MV , MW 
NODFIX( IFIX) -JNEW( NODF ) 

NU-JNEW( NODF ) *3-2 
NUNU-MSU1 1( NU )+MB( HU ) 

NVNV— MSUM( NU+-1)+MB(NU+1 ) 

NWNW-HSUM( NU+2 )+MB( NU+2 ) 

AA( NUNU) -AA( NUNU )+MU* SK2 
AA( NVNV ) - AA( NVNV )+MV*S K2 
AA( NWNW) -AA( NWNW)+MW* SK2 
MPTAB(NU)-MU 
MPTAB( NU+1 ) -MV 

416 UPTAB ( NU+2 ) - MW 

417 CONTINUE 
IF(NLOAD.LE.O) GOTO 739 
DO 41 IL-1, NLOAD 

READ(5 , 42) NODLOD(IL) ,PX,PY,PZ 
IZ-NODLOD(IL) 

WRITE(6 ,43) NODLOD(IL) ,PX,PY,PZ 
NODLOD(IL) -JNEW( IZ) 

IZ1— (JNEW(IZ)— 1)*3+1 
BB(IZ1 ,1)-PX 
BB( IZ1+1 , 1)-PY 
BB( IZ1+2 , 1)«PZ 

41 CONTINUE 

43 FQRMAT(5X, 15 ,3(F12 . 5 , 2X) ) 

739 • IF(NSPD.LE.O) GOTO 738 
DO 735 N— 1 , NSPD 
READ(5,736) NODS ,K,DISP(N) 

• WRITE(6,737) NODS,K,DISP(N) 

. NU-(JNEW(N0DS)-1)*3+K 
NDISP(N)-NU 
NUNU-MSUM( NU )+MB( NU ) 


AA(NUNU)«SK2 

735 BB(NU,1)-SK2*DISP(N) 

736 F0RMATC2I5.E14.5) 

737 / FORMATC5X,I5,2X,Il,3X,E12.5) 

738' CONTINUE 

R(l;NfiOF)-BB(l,l;NDOF) 

IZIP1-16 

CALL Q3CLOCKS( CPU , HALL) 

WRITE(6,9999) IZIPI, CPU, WALL 
IFAC-0 
ALP-0 
IZIPI-2X 

CALL SYMBAN( LNSTIF , NDOF , MB , MSUM , AA, 1 , BB , IF AC , T1 , IERR 
1,ALP,Z,T2,T3,T4,I) 

IZIP1-22 

CALL Q3CLOCKSC CPU, WALL) 

WRITE(6,9999) IZIPI, CPU, WALL 
IF( IERR.EQ. 1) WRITE(6 ,415) IERR 
415 F ORMAT ( / / 10X' IERR- ' 12 , 1 OX' NONPOS ITI VE DEFINITE MATRIX') 

IF(IERR.NE.O) STOP 

C PRINT OUT UNIT LOAD DISPLACEMENTS AND STRESSES. 

C 

9000 CONTINUE 

WRITE (6, 425) 

425 FORMAT( 1H1// / 10X, 'UNIT LOAD DISPLA4EMENTS AND STRESSES'//) 

WRITE (6, 41 8) 

418 F0RMAT(6X, 4HN0DE, 6X, 1HX, 15X, 1HY, 13X, 1HZ , 13X, 1HU, 13X, 

1 1HV.13X, 1HW/) 

DO 551 N-l ,HN0DE 
551 IIMAX(N)-(N-1)*3+1 

U ( 1 ; NNODE)- Q8VGATHRC BB( 1,1; NDOF) , IIMAX( 1 ; NNODE) ; U( 1 ; NNODE) ) 

I IMAX( 1 ; NNODE) -IIMAX( 1 ; NN0DE)+1 

V( 1 ; NNODE) - Q8VGATHR( BB( 1 , 1 ; NDOF) , IIMAX( 1 ; NNODE) ; V( 1 ;.NNODE ) ) 

IIMAX( 1 ; NNODE) -IIMAXC 1 ; NNODE )+l 

W(1 ; NNODE ) - Q8VGATHRC BB( 1 , 1 ; NDOF) , I IMAX( 1 ; NNODE ) j W( 1 ; NNODE ) ) 

DO 944 IN-1 .NNODE 

944 WRITE(6 , 420) IN.XR(IN) ,YR(IN) ,ZR(IN) ,U(IN) ,V(IN) ,W(IN) 

420 FORMAT(5X,I5,2X,3(2X,E11.5),3(2X,E12.4)) 

C 

c *** 

c 

PYLD-0.0 

WRITE(6,306) 

306 F0RMATC1H1//10X, 'ELASTIC STRESSES: SX, SY, SZ, AND SYZ, SZX, SXY') 

307 FORMAT (5X, 16, 2X, 6E12.4) 

IGAUSP-0 

DO 300 IE-1, NELM 
NCUBE(1 ;8)«M0DE(1 ,IE;8) 

DO 301 1-1,8 
Il-NCUBE(I) 

DISCI, l)-U(Il) 

DIS(I, 2)-V(Il) 

301 DISCI, 3)-WCIl) 

CALL CORDINC NCUBE , MXNOD , XR , YR, ZR, XE ) 

CALL STRESS C DIS , XE , D , STRV , STRSV , BMT , WDUM) 

IL0C-CIE-1)*NQD2SR+1 
XC IL0C;NQD2SR)«STRSVC 1 , 1;NQD2SR) 

Y C ILOC ; HQD2SR)— STRVC 1,1; NQD2SR) 

DO 350 IG-1.NQD2 
DO 360 J-1,6 

360 ' STRCJ)-STRSV(IG.J) 

‘ IGAUSP-IGAUSP+L 
• CALL SEQUCSTR.SEFF) 

SEFF-SEFF/SIGBARCIGAUSP) 

IF C PYLD . GT . SEFF ) GOTO 350 



PYLD-SEFF 
IEY-IE 
. IGAUSY-IG 
350 . CONTINUE 

WRIXE(6 ,307) IE,(STR(J) ,J-1,6) 

300 • CONTINUE 

PYLD-l./PYLD 

WRITE (6, 305) IEY,IGAUSY,PYLD 

305 F0RMAT(1H1/ /10X, 'ELEMENT#' ,15, 2X, 'GAUSS PT-' ,I2,2X,'L0AD FACTOR AT 
C 'YIELD-', E12.6) 

READ(5,450) NTYP, NLM, SCRIT, RP.ACURCY 
WRITE ( 6 , 4 I 2 ) NTYP , SCRIT , NLM , RP , ACURCY 
412 F0RMAT(//9X, 'CRACK GROWTH CRITERION NTYP-', 12,' AND CTOD -'.E10.4, 
C//10X, 'NUMBER OF INCREMENTS TO RELEASE CRACK TIP FORCE-', 12, 

C/ / 10X , ' RELAXATION PARAMETER- ' , F5 . 2 , ' ( NORMAL) ' , 

C//10X, 'CRACK OPENING DISPLACEMENT ACCURCY-' ,E12.4) 

450 F0RMAT(2I5 .4E10.4) 

IF(NEP.EQ.O) STOP 
CALL PLAS 
IZIP1-26 
9991 STOP 
END 

SUBROUTINE OPEN(LOPN) 

COMMON/MAIN/ AA(2000000) ,BB(9600,1) ,D(6,6) ,DINV(6,6) , 

1DISP(80) ,EPS(12400) ,EFEST(1550) ,F0RCE(10) ,LINE(100) , 

2L0CAT(10) ,LBF0R(10) ,MB(9600) ,MSUM(9600) ,MPTAB(9600) , 

3MPLAS( 12400) ,M0DE(8,1550) ,MPLC(1550) ,NODXO(70O) , 

4N0DY0(700) ,N0DZ0(700) ,NODXC(700) ,NODYC(700) ,N0DZC(700) , 

5N0DFIX(80) , NODL0D(80) ,NDISP(80) ,R(9600) , 

6SIGBAR( 12400) ,SK(24,24) ,T1(9600) ,12(9600) ,T4(2000) , 

7T3(9600) ,U(3200) ,UOLD(3200) ,V(3200) ,V0LD(3200) ,V2(3200) , 

8W(3200) , WOLD(3200) ,X(74400) ,XR(3200) ,Y(74400) , 

9YR(3200) ,Z(9600) ,ZR(3200) 

COMMON/ CNST/ EPS I , SK2 , LMAX , KMAX , DAX , 

1PYLD, SCRIT .YOUNG, POIS , CRACK, PT, WIDTH, PMAX, HP , 

2SBAR, LPRIT , NGAUS , NLAYER, NNODE , 

3IN0DX0 , INODYO , INODZO .INODXC , INODYC , INODZC .LNSTIF ,MXNOD , 

4MXNEL , MXGAUS , ICUT , LTOTB , ITNODX , KLU , NTYP , NLM , 

5ND0F .KNEW , NEP ,ERIT ,NELM, AM, ROM 
C 

C THIS SUB-PROGRAM OPENS ALL THE Q30PNMAP FILES 

C MAXIMUM LENGTH OF ANY FILE IS 5376 SMALL PAGES (DECIMAL) 

C 

C TO CHANGE THE MAXIMUM LENGTH CHANGE THE DATA CARD 

C DATA LMAX / / 

C 

CHARACTER* 8 FILE, U0RD(8) 

DATA WORD/ 'ASTIF001' , 'ASTIF002' , 

Z 'ASTIF003' , 'ASTIF004' , 

Z 'ASTIF005' , 'ASTIF006' , 

Z 'ASTIF007' , 'ASTIF008' / 

DATA LMAX / 5376/ 

IF(LOPN.LE. LMAX) GO TO 20 
LOPNA- LMAX 
LDUM- LOPN/LMAX 
LOPNB- L0PN-LDUM*LMAX 
DO 10 I-l.LDUM 
FILE- WORD(I) 

ISTART- LMAX*512*(I-1)+1 

• CALL Q30PNMAP (IERR, FILE, AA( ISTART) , LOPNA, 1) 

PRINT 100, IERR, FILE, LOPNA 
WRITE(6, 100) IERR, FILE, LOPNA 
‘ IF(IERR.NE.O) STOP 

100 FORMAT (10X,' IERR FROM 0PN11AP-' ,Z16 ,5X, ' FILE ',A8, 

Z 2X,' IS OF LENGTH ',I10,2X,' SMALL PAGES (DECIMAL)',/) 

10 CONTINUE 


IF( LOPNB. EQ.O) RETURN 
ISTARX- LMAX*LDUM*5 12+1 
FILE- WORD( LDUM+1) 

‘ CALL Q30PNHAP (IERR, FILE, AA( I START) , LOPNB, 1) 

' PRINT 100, IERR, FILE, LOPNB 
WRITE ('6 , 100) IERR, FILE, LOPNB 
IF(IERR.NE.O) STOP 
RETURN 

20 CONTINUE 

FILE- WORD(l) 

CALL Q30PNMAP ( IERR, FILE, AA(1), LOPN, 1) 

PRINT 100, IERR, FILE, LOPN 

WRITE(6, 100) IERR, FILE, LOPN 

IF(IERR.NE.O) STOP 

RETURN 

END 

FUNCTION FNMAT(SBAR, CE, EPS ,M) 

COMMON/ MLTNMAT/YSTRS( 20) ,YSTRN(20) ,PLMODR(20) ,NSEGMT 
C 

C EPST- TOTAL STRAIN 

C EPS - PLASTIC STRAIN 

C 

EPST-EPS+SBAR/CE 
DO 10 I-l.NSEGMT 

10 IF(EPST.LT. YSTRN(I) ) GOTO 11 

11 FNMAT»PLMODR(I)*CE 
RETURN 

END 

SUBROUTINE HMPLN ( AA.HSUM , MB , MPTAB , NODP , INOD , SK2 , ID , NDOF , LNSTIF ) 
C 

C *** IMPOSING SYMMETRIC BOUNDARY CONDITIONS 
C 

DIMENSION AA( LNSTIF) ,MSUH(NDOF) ,MB( NDOF) .MPTAB(NDOF) .NODP(INOD) 

DO 100 I- 1, INOD 

L-NODP(I) 

NU-(L-1)*3+ID 
NUNU»MSUM( NU)+MB( NU) 

MPTAB(NU)-NU 

100 AA(NUNU)-AA(NUNU)+SK2 
RETURN 
END 

SUBROUTINE NSAMC(MSAME , NSAUE, IIMAX,MB,LHAX,KMAX,NODPEL,MXNOD, 
1MXNEL.ND0F) 

DIMENSION MSAME(NODPEL.MXNEL) ,N3AME(MXN0D,20) .IIMAX(MXNOD) , 

C MB(NDOF) 

£ ******* A***************************************** ********* 

C MXNEL - MAXIMUM NUMBER OF ELEMENTS 

C MXNOD - MAXIMUM NUMBER OF NODES 

C NODPEL - It OF NOED PER ELEMENTS 

C LMAX - It OF NOEDS IN THE PROBLEM 

C KMAX - it OF ELEMENTS IN THE PROBLEM 

C MSAME(NODPEL, IEL) - NODEL CONNECTIVITY IF I EL ELEMENT 
C NDOF - LMAX* it OF DOF PER NODE 

C MB(NDOF) » BAND WIDTHS OF ALL NDOF DEGREE-OF- FREEDOM 

C 

c ********************************************************** 

DO 10 IE-1, KMAX 
DO 20 J-l, NODPEL 
IK-MSAME(J.IE) 

. IIMAX(IK)-IIMAX(IK)+1 
20 NSAME(IK,IIMAX(IK) )-IE 

10 ' CONTINUE 

C16 ' FORMAT(16I5) 

C . ***CALCULATE MB VECTOR 

IBANDW-0 

DO 350 NODE-1 ,LHAX 
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MAXDIF-0 

IM-IIMAX(NODE) 

. DO 351 M-l.IM 
. NTRI-NS AME ( NODE , M) 

DO 351 L-l.NODPEL 
NUM-MSAME(L,NTRI) 

NDIFF-3* ( HUM-NODE ) 

IF(NDIFF.LT.IIAXDIF) MAXDIF-NDIFF 
351 CONTINUE 

IF( IBANDW.LT. IABS(MAXDIF) ) IBANDW-IABS(HAXDIF) 

NU-3* ( NODE- 1 )+l 

NV-NU+1 

NW-NV+1 

MB(NU)»IABS(MAXDIF)+1 
IF(MB(NU) .GT.NU) MB(NU)-NU 
MB(NV)-MB(NU)+1 
MB(NW)-MB(NV)+1 
350 CONTINUE 

IBANDU-IBANDW+3 
WRITE (6 ,25) IBANDW 

25, FORMAT (5X, 'MAX BAND WIDTH- ' ,16) 

RETURN 

END 

SUBROUTINE DC0N( YOUNG , POIS , D , DINV) 

*** 3_d D(6,6) & DINV MATRICES FOR ISOTROPIC MATERIAL 

DIMENSION D(6 ,6) ,DINV(6 ,6) 

DEL-YOUNG* ( 1-POIS) / ( ( 1+POIS) * ( 1-2*P0IS) ) 
DEL2-P0IS/ (1-POIS) 

DEL3“(1-2*P0IS)/ (2*(1-P0IS) ) 

D(l,l;36)-0.0 
DINV (1,1.; 36) -0.0 
D(1,1)-DEL 
D(1,2)-DEL*DEL2 
D( 1 ,3)-D(l , 2) 

D(2 , 2)-D(l , 1) 

D(2 ,3)-D(l ,3) 

D(3,3)-D(l,l) 

D(4,4)-DEL*DEL3 

D(5,5)-D(4,4) 

D(6,6)-D(5,5) 

*** INVERSE OF D-HATRIX 
DINV(1,1)-1. /YOUNG 
DINV (1, 2)— POIS/ YOUNG 
DINV(1 ,3)-DINV(l ,2) 

DINV(2,2)-DINV(1,1) 

DINV(2 , 3)-DINV( 1 ,2) 

DINV(3,3)-DINV(1,1) 

DINV(4,4)“ 2. *( 1+POIS) /YOUNG 

DINV(5,5)-DINV(.4,4) 

DINV(6 ,6)-DINV(4 ,4) 

DO 5 1-1,3 
DO 5 J-I,3 
D(J,I)-D(I,J) 

DINV(J,I)-DINV(I,J) 

CONTINUE 

RETURN 

END 

SUBROUTINE SHAPE(X,Y,Z,R) 

*** SHAPE FUNCTIONS 

DIMENSION R(8) 

R(l)-1. 

R(2)-X 
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R(3)-Y 

R(4)-Z 

R(5)-X*Y 

R(6)«Y*Z 

R(7)-Z*X 

• R(8)-X*Y*Z 
RETURN 
END 

SUBROUTINE CORDIN( NCUBE , MXNOD, XR , YR , ZR , A) 

*** EVALUATE A(8,3) CARTESIAN COORDINATE MATRIX 

DIMENSION A(8 ,3) ,NCUBE(8) ,XR(MXNOD) ,YR(MXNOD) ,ZR( MXNOD) 

DO 1 1-1,8 
Nl-NCUBE(I) 

A(I, l)-XR(Nl) 

A(I,2)-YR(N1) 

1 A( I, 3)-ZR(Nl) 

RETURN 

END 

SUBROUTINE ACAL 
COMMON/ AINV/AI(8, 8) 

COMMON/GENRL/GCR(8 ,3) 

DIMENSION Rl(8) ,DUM(8,1) ,IPIVOT(8) ,IHK(16) ,A2(8,8) 
A2(l,l;64)-0.0 
DO 1 1-1,8 
Xl-GCR(I.l) 

Y1-GCR(I,2) 

Z1-GCR(I,3) 

CALL SHAPE(X1,Y1,Z1,R1) 

DO 1 J-1,8 

1 A2(I, J)’-Rl(J) 

CALL MATINV(A2,8,8,DUM,1,0,DET) 

AI,(1,1;64)-A2(1,1;64) 

RETURN 

END 

SUBROUTINE DERIVE (X,Y,Z,R) 

COMMON/ AINV/AI(8, 8) 

ROWWISE DN(3 , 8) 

DIMENSION R(3,8) 

DN(1 , 1; 24)-0 . 0 
C**** DN/DXI NOW 

DN( 1 ,2)-l . 0 
DN(1,5)-Y 
DN(1,7)-Z 
DN(1,8)-Y*Z 

C**** dn/deta now 

DN(2,3)“1.0 
DN(2,5)-X 
DN(2,6)-Z 
DN(2 ,8)- X*Z 
C**** DN/DZETA NOW 

DN(3 ,4)-l .0 
DN(3,6)-Y 
DN(3,7)-X 
DN(3 ,8)-X*Y 
DO 10 J-1,3 
DO 10 1-1,8 

R(J,I)- Q8SD0T ( DN(J,1;8) , AI(l,Ij8) ) 

10 . CONTINUE 

RETURN 
' END 

• SUBROUTINE SMALLK( SMK, XE, 'D, IERR) 

C 

C THIS MODULE GENERATES AN ELEMENTAL STIFFNESS MATRIX FOR THE GIVEN 
C ELEMENT. VECTOR VERSION 
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DIMENSION SMK(24,24), XE(8 ,3) ,D(6 ,6) 

, COMMON/ D382/NNPE,HDF,NQD,NSTR,NQD2,NNPE2,NQD2NPE,NQD2SR,MXQ2S 
, READ KE 

DIMENSION KE(324) 

C • ' 

C 

C INPUT: D(6,6) - MODULUS MATRIX 

C XE(8,3) - 8 NODES X,Y, Z COORDINATES 

C OUTPUT : SMK(24,24) - STIFFNESS MATRIX 

C 

C 

C 

C 

C KE - ELEMENTAL STIFFNESS MATRICES FOR ALL DISTINCT ELEMENTS, IN 

C ROWISE NODAL BLOCK LOWER TRIANGULAR FORM 

C PE - ELEMENTAL LOAD VECTORS FOR ALL ELEMENTS IN NODAL BLOCK FORM 

DATA NDFX, NSTRX, NNPEX / 3, 6, 8 / 

C 

C 

C -i 

C 

C NDF - NUMBER OF DISPL. DEGREES OF FREEDOM PER NODE 

C NSTR - NUMBER OF STRESS RESULTANTS PER NODE 

C NQD - NUMBER OF QUADRATURE POINTS IN EACH DIRECTION 

C NNPE - NUMBER OF NODES PER ELEMENT 

C 
C 

C*** ETH- THERMAL STRAINS IN THE CARTESIAN SYSTEM. 

C**** FTHERM- THERMAL LOAD VECTOR. 


C 

C 

C [D] - STRESS STRAIN MATRIX 
C 

DIMENSION IBSP(6 , 3) , B(64,3), BJ(288,3) ,CK(288,3) , 

Z WTDETEX(64) , SUM(288) 

DIMENSION UTDET(8) 

DIMENSION CTH(64) TPST(8) 

DATA IBSP/ 1, 2*6, >, 0, 3, 

Z 0, 2, 0, i, 3, 0, 

Z 2*0, 3, 0, 2, 1 / 

C 

C [IBSP] - SPARSITY PATTERN AND POINTER MATRIX FOR [B] AND [BJ] 

C [B] - STRAIN DISPLACEMENT MATRIX 

C [BJ] - ANOTHER STRAIN DISPLACEMENT MATRIX 

C [CK] - A ROW FOR EACH STIFFNESS MATRIX NODAL PARTITION 

C (WTDETEX) - REPLICATED S/EIGHTED DETERMINANTS 

C (SUM) - TEMPORARY STORAGE 

C 


DIMENSION IREPL(36) ,IP0SN(210) 

DESCRIPTOR IREPLD, SORCD, DESTD 
C 

C IREPL 
C 

C IPOSN 
C 

C IREPLD 
C SORCD 
C DESTD 
C 

DATA LENI, LENB, LENBJ, LENC, LENW, LENWT 
Z / 18, 192, 864, 864, 288, 64 / 

C 

C THESE ARE THE DIMENSIONED LENGTHS OF [IBSP], [B], [BJ], [CK] , 
C (WTDETEX), AND (SUM) FOR ZEROING OUT PURPOSES. 

C 


- VECTOR OF LENGTH "NNPE2” CONTAINING ZEROS USED IN THE 
REPLICATION PROCESS 

- ARRAY OF LENGTH "NNPE2" USED TO CORRECTLY POSITION 
THE NODAL PARTITIONS IN [KE] 

- VECTOR DESCRIPTOR FOR (IREPL) 

- VECTOR DESCRIPTOR FOR THE REPLICATION SOURCE 

- VECTOR DESCRIPTOR FOR THE REPLICATION DESTINATION 


l 


i'' 



i 


h 
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L3 - NQD2NPE - NQD2 + 1 
DO 140 KK - 1, NNPE 

ASSIGN IREPLD, IREPL(1; KK) 

ASSIGN SORCD, B(L3,IT; NQD2) t 

ASSIGN DESTD, BJ(L1,IT; L2) 

. • * CALL Q8VXTOV(X'02' , 0, IREPLD, 0, SORCD, 0, DESTD) 

L2 » L2 + NQD2 

LI - LI - L2 

L3 “ L3 - NQD2 

140 CONTINUE 

150 CONTINUE 

C 

DO 155 IT-1,3 

B(1,IT;NQD2NPE) - B ( 1 , IT ; NQD2NPE ) *WTDETEX ( 1 ; NQD2NPE ) 
155 CONTINUE 


FORM ( [ B ] **T * [D]) * [BJ] A ROW AT A TIME. 


L9-1 

DO 400 II - 1, NDF 
CK(1,1; LENDF) - 0. 

DO 230 KK - 1, NSTR 
SUM(1; NQD2NPE) - 0. 

DO 210 JJ - 1, NSTR 
IT - IBSP(JJ,II) 

IF (IT .EQ. 0) GOTO 210 
IF (D(JJ,KK) .EQ. 0.) GOTO 210 

SUM(1; NQD2NPE) - SUM(1; NQD2NPE) + B(1,IT; NQD2NPE) 

* D(JJ,KK) 


210 CONTINUE 


FILL UP THE REST OF (SUM) 


IF (SUM(l) .EQ. 0.) GOTO 225 
LI - LE - NQD2 + 1 
L2 - NQD2 

L3 - NQD2NPE - NQD2 + 1 
DO 215 JJ - 2, NNPE 

SUM(L1; L2) - SUM(L3; L2) 

L2 - L2 + NQD2 
LI - LI - L2 
L3 - L3 - NQD2 
215 CONTINUE 

DO 220 JJ - 1, NDF 
IT - IBSP(KK.JJ) 

IF (IT .EQ. 0) GOTO 220 

CK(1,JJ ; LE) - CK(1,JJ ; LE) + SUtl(l; LE) * BJ(1,IT; LE) 
220 CONTINUE 

225 CONTINUE 

230 CONTINUE 

C WE NOW HAVE THE II-TH ROW (BEFORE SUMMING) FOR ALL 
C "NNPE2" NODAL PARTITIONS OF THE ELEMENTAL STIFFNESS MATRIX. 

C 

DO 310 JJ - 1, NDF 
LI - 1 

DO 300 KK - 1, NNPE2 
L2 - L9 + IPOSN(KK) 

KE(L2) - Q8SSUM(CK(L1,JJ; NQD2)) 

LI - LI + NQD2 
300 CONTINUE 

L9 - L9 + 1 
310 • CONTINUE 
400 . CONTINUE 
C 

L9-1 
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Z , (B(6401) ,DAJ23(1)) 

Z , (B(7681) ,DAJ31(1)) 

Z , ( B(8961) ,DAJ32(1) ) 

Z , (B(1024l) ,DAJ33(1)) 

Z ,(B('11521) ,AJU (1 )) 

. Z -,(8(11585), AJ12 (1 )) 

Z ,(B(11649),AJ13 (1 )) 

Z ,(B(11713) , AJ21 (1 )) 

Z ,(B(11777) ,AJ22 (1 )) 

Z ,(B(11841),AJ23 (1 )) 

Z ,(B(11905) ,AJ31 (1 )) 

Z ,(B(U969),AJ32 (1 )) 

Z , (B( 12033) ,AJ33 (1 )) 

Z ,(8(12097) ,AJI11(1)) 

EQUIVALENCE (B( 12161) ,AJI12(1 ) ) 

Z ,(B(12225) ,AJI13(1)) 

Z ,(B( 12289) ,AJI21(1)) 

Z ,(B(12353) ,AJI22(1)) 

Z ,(B(12417) ,AJI23(1)) 

• Z ,(B( 12481) ,AJI31(1) ) 

’ Z ,(B(12545) ,AJI32(1) ) 

Z , (B( 12609) , AJI33(1) ) 

Z ,( B( 12673) ,DET11( 1) ) 

Z ,(B(12737),DET12(1)) 

Z , ( B( 12801) ,DET13(1) ) 

Z ,(B(12865) ,DET21(1)) 

Z ,(B( 12929) ,DET22(1)) 

Z , (B( 12993) ,DET23(1) ) 

Z ,(B( 13057) ,DET31(1) ) 

Z ,(B(13121) ,DET32(1)) 

Z ,( B(13185) ,DET33(1) ) 

C DIMENSION BJ(216,3) 

C**** THE ABOVE DIMENSIONS ALLOW UPTO 4 POINT GAUSSIAN IN EACH DIRECTION 
LEN-11520 , 

LENI-13248 

CALL ZEROLV ( BB , LEN ) 

CALL ZEROLV ( B , LENI ) 

NG- N*N*N 

C*** NG ARE THE TOTAL NO OF INTEGRATION POINTS 
NS-8 

C**** NS- NO OF SHAPE FUNCTIONS 
NSTR-6 
NCORD-3 
NFREE-3 

C*** NSTR- NO OF STRAINS. NCORD- NO OF COORDINATES NFREE- NO OF DOF P 
MAX-NG*NSTR 
DO 10 I-l.N 

X- CORD( I,N) 

XI- (X+l.)/2. 

WI-WEIGHT(I.N) 

DO 10 J»1,N 
Y- CORD(J.N) 

ETA-(Y+l.)/2. 

UJ- WEIGHT(J,N) 

DO 10 K -1,N 
Z- CORD(K,N) 

ZI-(Z+1.0)/2.0 

WK-UEIGHT(K.N) 

CALL DERIVE ( XI, ETA, ZI, R) 

• II-N*N*(I-1)+N*(J-1)+K 
W(II)-WI*WJ*WK/8.0 

DO 20 IJ-l.NS 
' IN-NS*(II-1)+IJ 

• DNX(IN)-R(1,IJ) 

DNE( IN)-R(2 , IJ) 

DNZ(IN) -R(3,IJ) 


20 CONTINUE 
10 CONTINUE 

C***** NOW GENERATE THE MASTER VECTOR OF THE COORDINATES 
, DO 30 J-l.NG 
NT-NS*(J-1)+1 
• XX(NT;tfS)-XE(l,l;NS) 

YY(NT;NS)»XE(1,2;NS) 

ZZ(NT;NS) - XE(1,3;NS) 

30 CONTINUE 

LN-NG*NS 

DAJ11(1;LN)-DNX(1;LN)*XX(1;LN) 

DAJ12( 1 j LN) -DNX( 1 ; LN)*YY( 1 ; LN) 

DAJ13(1;LN) - DNX(1;LN)* ZZ(1;LN) 

DAJ21(1;LN)«DNE(1;LN)*XX(1;LN) 

DAJ22(1;LN)-DNE(1;LN)*YY(1;LN) 

DAJ23( 1 ;LN) - DNE(1;LN)* ZZ(L;LN) 

DAJ31(1;LH) - DNZ(ljLN)* XX(1;LN) 

DAJ32(1;LN) - DNZ(1;LN)* YY ( 1 ; LN ) 

DAJ33(1;LN) - DNZ(1;LN)* ZZ(ljLN) 

DO 40 I-l.NG 
NT-NS*(I-1)+1 

AJll(I)- Q8SSUM(DAJ11(NT;NS) ) 

Q8SSUl'l(DAJ12(NT;NS) ) 

«Q8SSUM(DAJ13(NT;NS)) 

Q8SSUM(DAJ21(NT;NS)) 

Q8SS UM( DAJ2 2 ( NT ; NS ) ) 

“Q8SSUM( DAJ23(NTj NS) ) 

“Q8SSUM(DAJ31(NT; NS) ) 

-Q8SSUM( DAJ32 ( NT ; NS ) ) 

-Q8SSUM(DAJ33(NT;NS)) 


AJ12(I)- 
AJ13(I) 

AJ2KI)- 
AJ22(I)- 
AJ23(I) 

AJ31(I) 

AJ32(I) 

AJ33(I) 

40 CONTINUE 

DET11(1;NG) “AJ22(1 ;NG)*AJ33( 1;NG)-AJ32(1 ;NG)*AJ23(1;NG) 

DET12(1 ;NG) “AJ21(1;NG)*AJ33(1;NG)-AJ31(1 ;NG)*AJ23(1 ;NG) 

DET13(1 ;NG) ”AJ21(1 ;NG)*AJ32(1 ;NG)-AJ31(1;NG)*AJ22(1 ;NG) 
DET21(1;NG) -AJ12(1;NG)*AJ33(1;NG)-AJ32(1;NG)*AJ13(1;NG) 

DET22( 1 ;NG) -AJ11(1;NG)*AJ33(1;NG)-AJ31(1;NG)*AJ13(1;NG) 

DET23(1 ;NG) »AJU(1;NG)*AJ32(1;NG)-AJ31(1;NG)*AJ12(1;NG) 
DET31(1;NG) -Aji.2(l;NG)*AJ23(l;NG)-AJ22(l;NG)*AJ13(l;NG) 
DET32(1;NG) -AJU( 1 ;NG)*AJ23( 1 ;NG)-AJ21( 1 ;NG)*AJ13(1 ;NG) 
DET33(1;NG) -AJU(1;NG)*AJ22(1;NG)-AJ21(1;NG)*AJ12(1;NG) 

DET(1;NG) «AJU(1;NG)*DETU(1;NG)-AJ12(1;NG)*DET12(1;NG)+ 
ZAJ13(1;NG)* DET13(1;NG) 

AJI11(1; NG) - DET1 1 ( 1 ; NG) / DET ( 1 ; NG ) 

AJI12(1; NG) — DET21(1;NG)/DET(1;NG) 

AJI13(1; NG) - DET3 1(1; NG) / DET ( 1 ; NG) 

AJT21(1; NG) — DET12 ( 1 ;NG) /DET(1 ;NG) 

AJI22(1; NG) - DET22(1;NG)/DET(1 ;NG) 

AJ 123(1; NG) — DET32( 1 ; NG) /DET ( 1 ; NG) 

AJI31(1; NG) - DET13(1;NG)/D£T(1;NG) 

AJI32(1; NG) — DET23(1;NG)/DET(1;NG) 

AJ 133(1; NG) " DET33(1;NG)/ DET(1;NG) 

C*** JOCOBIANS AND THEIR INVERSES ARE READY 
DO 50 J»1,NG 
NT«NS*(J-1)+1 

DSX(NT;NS) -DNX(NT;NS)*AJI11(J)+DNE(NT;NS)*AJI12(J) +DNZ(NT;NS)* 
Z AJI13(J) 

DSY(NT;NS) 

Z AJI23( J) 

DSZ(NT;NS) 

• Z AJI33(J) 

50 CONTINUE 

C***'* CARTESIAN DERIVATIVES ARE READY 
' IF(IC0DE.EQ.2) GO TO 320 
- 11-1 
I2-NS 

CALL Q8INTVAL (0,0, II ,0, 12 ,0,INVA( 1 ;NG) ) 


-DNX(NT;NS)*AJI21(J) +DNE(NT;NS)*AJI22(J) +DNZ(NT;NS) 
”DNX(NT;NS)*AJI31( J) +DNE(NT;NS)*AJI32(J) +DNZ(NT;NS) 



LN- NS*NG 
DO 60 I-l.NS 
NT- NG*(I-l)+l 

■ DNX( NT ; NG) -Q8V GATIIR( DSX( X ; LN ) ,INVA(1 ;NG) ;DNX(NT;NG) ) 

DNE '( NT ; NG) -Q8 VGAT1IR( DSY ( I ; LN) , INVA( 1 ; NG ) ; DNE ( NT ; NG) ) 

. DNZ(NT-NG) -Q8VGATHR(DSZ(I;LN) , INVA( 1 ;NG) ; DNZ(NT;NG) ) 

60 CONTINUE 

C**** CARTESIAN DERIVATIVES ARE REORDERED SO THAT THE DERIVATIVES AT A 
C*** GAUSSIAN POINTS ARE GROUPED 
LEN-NS*NG 

BJ(1,1;LEN)“ DNX(1 ;LEN) 

BJ(L,2;LEN)- DNE(1; LEN) 

BJ(1,3;LEN)“ DNZ(ljLEN) 

C***** COMPUTE THE PRODUCT OF WEIGHT AND DETRMINENTS 

WDUM( 1 ;NG)- DET(1 ;NG)*W(1 ;NG) 

RETURN 

320 CONTINUE 

LEN“NS*NG 

BJ(1,1;LEN)- DSX(1;LEN) 

BJ(1,2;LEN)“ DSY(1;LEN) 

BJ(1 , 3;LEN)- DSZ(I;LEN) 

WDUM(1;NG)— DET(l;NG)*W(l;NG) 

RETURN 

END 

SUBROUTINE STRESS(DIS ,XE , D , STR, SIRS , B , WDUM) 

C 

C 

COfIMON/D382/NNPE,NDF,NQD,NSTR,NQD2,NNPE2,NQD2NPE,NQD2SR,MXQ2S 
DIMENSION DIS(8,3), XE(8,3) ,D(6 , 6) 

DIMENSION STRS(8 , 6) ,STR(8 ,6) 

DIMENSION HDUM( 8), SUM(64) , B( 64,3), DISP(64,3) 

DIMENSION IREPLC20) ,STRD(6,'8) 

DESCRIPTOR - IREPLD , SORCD , DESTD 
C 

DIMENSION IBSP(6 ,3) 

DATA IBSP / 1,2*0, 2, 0, 3, 

Z 0, 2, 0, 1, 3, 0, 

Z 2*0, 3, 0, 2, 1 / 

C 

c 

DATA NS, NSTR, NSH, NFREE / 8, 6, 3, 3 / 

c ******* NS- NUMBER OF SHAPE FUNCTIONS OR NODES ON THE ELEMENT 
c ******* NSTR- NUMBER OF STRAINS 

c******* NSH— NUMBER OF INDEPENDENT DERIVATIVES IN THE B MATRIX 

C******* NFREE- NUMBER OF DEGREES OF FREEDOM PER NODE 

C 
C 

LEN- NSTR*NQD2 
LDISP-NQD2NPE*NSH 
IREPL(l;NQD2)-0 
STRS( I , 1; LEN) -0.0 
CALL ZEROLV( DISP , LDISP) 

C 

C**** PICK UP THE U, V,W DISPLACEMENTS SEPARATELY 
C 

C**** REPLICATE THE DISPLACEMENTS NQD2 TIMES 
C 

ASSIGN IREPLD, IREPL(1;NQD2) 

DO 25 KC-1, NFREE 
• ASSIGN SORCD, DIS(1,KC;NS) 

ASSIGN DESTD , DISP(1,KC;NQD2NPE) 

CALL Q8VXT0V (X'02' , 0, IREPLD, 0, SORCD, 0, DESTD) 

25 ‘ CONTINUE 

c **** THE MASTER DISP VECTOR READY 
C 

C**** GET THE CARTESIAN DERIVATIVES AT THE NODES 
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NJ-(J-1)*HQD2+1 

220 F0RC(IJ)-Q8SSUM(SUM(NJ;NQD2)) 

400 CONTINUE 
‘ RETURN 
END 

. • SUBROUTINE MATINV(A,NMAX, N,B, MAX, M, DETERM) 
DIMENSION A(NMAX,NMAX) ,B(NMAX,MAX) 

DIMENSION IPIVOT(IOO) ,INDEX(100,,2> ,PIV0T(100) 

IF M-0 IT CALCULATES THE INVERSE ONLY. 

IF M-l IT CALCULATES THE SOL TO AX-B IN B 
INITIALIZATION 

10 DETERM-1.0 
15 DO 20 J-l ,N 
20 IPIV0T(J)-0 
30 DO 550 1-1 ,N 

SEARCH FOR THE PIVOT ELEMENT 

40 AMAX-0.0 
45 DO 105 J-1,N 
50 IF(IPIVOT(J)-1)60,105 ,60 

60 DO 100 K-l ,N 
70 IF(IPIV0T(K)-1) 80, 100,740 
80 IF( ABS(AMAX)- ABS( A( J,K) ) )85 , 100 , 100 
85 IROW- J 
90 ICOLUM-K 
95 AMAX- A(J,K) 

100 CONTINUE 
105 CONTINUE 

110 IPIV0T(IC0LU11)-IPIV0T(IC0LUM)+1 

INTERCHANGE ROWS TO PUT ELELMENG ON DIAGONAL 

130 IF( IRGW-ICOLUM) 140 , 260 , 140 
140 DETERM- -DETERM 
150 DO 200 L— 1 ,N 
160 SWAP- A( IROW.L) 

170 A(IR0W,L)-A(IC0LUM,L) 

200 A(IC0LUM,L)- SWAP 
205 IF(M)260,260,210 
210 DO 250 L-l.M 
220 SWAP- B( IROW.L) 

230 B( IROW.L)- B( ICOLUM, L) 

250 B(IC0LUM,L)- SWAP 

260 IND£X( I, 1)» IROW 

270 INDEX(I , 2)» ICOLUM 

310 PIVOT(I)- A(IC0LUM, ICOLUM) 

120 DETERM- DETERM* PIVOT( I) 

DIVIDE 'PIVOT BY PIVOT ELEMENT 

330 A( ICOLUM, ICOLUM) -1.0 
340 DO 350 L-l ,N 

350 A(ICOLUM,L)“ A(ICOLUM,L)/PIVOT(I) 

355 IF(M) 380,380,360 
360 DO 370 L-l ,H 

370 B(IC0LUtl,L)« B(ICOLUM,L)/PIVOT(I) 

REDUCE NON-PIVOT ROWS 

380 DO 550L1-1 ,N 

390 IF(Ll- f ICOLUll)400 , 550 , 400 

400 T- A(L1, ICOLUM) 

420 A( LI , IC0LUlI)-0 . 0 
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LN- NS*NG 
DO 60 I-l.NS 
NT- NG*(I-1)+1 

DNX(NT ; NG)-Q6VGATHR(DSX( I ; LN) , INVA( 1 ; HG) ; DNX( NT; NG) ) 

DNE(NT ; NG) -Q8VGATUR(DSY ( I ; LN) ,INVA(1.;NG) ;DNE(NT;NG) ) 

DNZ(NT-NG) -Q8VGATHR(DSZ(I;LN),INVA(1;NG);DNZ<NT;NG)) 

60 CONTINUE 

C**** CARTESIAN DERIVATIVES ARE REORDERED SO THAT THE DERIVATIVES AT A 
C*** GAUSSIAN POINTS ARE GROUPED 
LEN-NS*NG 

BJ(I,1;LEN)» DNX(1;LEN) 

BJ(1,2;LEN)« DNE( 1; LEN) 

BJ(1,3;LEN)- DNZ(1;LEN) 

C***** COMPUTE THE PRODUCT OF WEIGHT AND DETRMINENTS 
WDUM( I ;NG)- DET(1;NG)*W(1;NG) 

RETURN 

320 CONTINUE 

L£N-NS*NG 

BJ(1,1;LEN)“ DSX(I;LEN) 

BJ( 1 , 2;LEH)« DSY(1;LEN) 

BJ(1,3;LEN)« DSZ(1;LEN) 

WDUll(ljNG)- DET(1;NG)*W(1;NG) 

RETURN 

END 

SUBROUTINE STRESS(DIS,XE,D,STR, STRS ,B,WDUM) 

C 

C 

C0IIM0N/D382/NNPE , NDF , NQD , NSTR, NQD2 , NHPE2 , NQD2NPE , NQD2SR, MXQ2S 
DIMENSION DIS(8,3), XE(8,3) ,D(6,6) 

DIMENSION STRS(8 , 6) ,STR(8 , 6) 

DIMENSION WDUM( 8), SUM(64) , B( 64,3), DISP(64,3) 

DIMENSION IREPL(20) ,STRD(6,8) 

DESCRIPTOR I RE PLD, SORCD , DESTD 
C 

DIMENSION IBSP(6 , 3) 

DATA IBSP / 1,2*0, 2, 0, 3, 

Z 0,2, 0,1, 3,0, 

Z 2*0, 3, 0, 2, 1 / 

c 
c 

DATA NS, NSTR, NSH, NFREE / 8, 6, 3, 3 / 

c ******* NS- NUMBER OF SHAPE FUNCTIONS OR NODES ON THE ELEMENT 
C******* NSTR- HUMBER OF STRAINS 

c******* NSH" NUMBER OF INDEPENDENT DERIVATIVES IN THE B MATRIX 

c ******* NFREE- NUMBER OF DEGREES OF FREEDOM PER NODE 

C 
C 

LEN- NSTR*N(JD2 

LDISP-NQD2NPE*NSH 

IREPL(1;NQD2)«0 

SIRS ( 1 , l;LEN)-0.0 

CALL ZEROLV(DISP.LDISP) 

C 

c **** PICK UP THE U,V,W DISPLACEMENTS SEPARATELY 
C 

c **** REPLICATE THE DISPLACEMENTS NQD2 TIMES 
C 

ASSIGN IREPLD, IREPL( 1 ;NQD2) 

DO 25 KC-I, NFREE 
• ASSIGN SORCD, DIS(1,KC;NS) 

ASSIGN DESTD , DISP( 1 ,KC;NQD2HPE) 

CALL Q8VXT0V (X'02', 0, IREPLD, 0, SORCD, 0, DESTD) 

25 ' CONTINUE 

C**** THE MASTER DISP VECTOR READY 
C 

c **** G E t THE CARTESIAN DERIVATIVES AT THE NODES 
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CALL CDER( XE,NyD, B, UDUM, 2) 


C****** NOW DO THE PRODUCT D * B* DISPLACEMENTS 
. DO 100 I»1,NSTR 
sum(1;Nqd2npe)-o.o 

■ DO 110'j-l,NSH 
XT- IBSP(I,J) 

IF(IT.EQ.O) GO TO 110 

SUM(1;NQD2NPE)« SUM(1;NQD2HPE)+B(1,IT;NQD2NPE)* DISP(1,J;NQL PE) 
110 CONTINUE 

DO 120 J-1.NQD2 
I1-(J~1)*NS+1 

STR(J,I)- Q8SSUM( SUH( II ; NS) ) 

STRD(I,J)-STR(J,I) 

120 CONTINUE 
100 CONTINUE 
C 

DO 130 I-1.NQD2 
DO 140 J-l.KSTR 

STRS(I, J)- Q8SD0T (D(l , J;NSTR) , STRD(1 ,1 ;NSTR) ) 

140 CONTINUE 
130 CONTINUE 
RETURN 
END 

SUBROUTINE FORCEP( BB , UTDET , STRS , FORC) 

C0M110N/D382/NNPE,NDF,NQD,NSTR,NQD2,NNPE2,NQD2NPE,NQD2SR,MXQ2 
DIMENSION IBSP(6,3) ,B(64,3) ,WTDET(8) ,STRS(8,6) 

DIMENSION WTDETEX(64) ,SUM(64) ,SIG(64,6) ,IREPL(20) 

DIMENSION FORC( 24) , BB(64 ,3) ,INDX(8) ,SH(8) 

DATA IBSP / 1,2*0, 2, 0, 3, 

Z 0, 2, 0, 1, 3, 0, 

Z 2*0, 3, 0, 2, 1 / 

C 

C NQD2-NQD*NQD*NQD 

C NQD2NPE-NQD2*NNPE 

DESCRIPTOR IREPLD, SORCD, DESTD 
DESCRIPTOR BDESC 
.IREPL(l;NNPE)-0 
ASSIGN IREPLD, IREPL(1; NNPE) 

ASSIGN SORCD, WTDET( 1 J NQD2) 

ASSIGN DESTD, WTDETEX(1; NQD2NPE) 

CALL Q8VXT0V(X'02', 0, IREPLD, 0, SORCD, 0, DESTD) 

DO 155 IT-1, NDF 
DO 156 J-l.NNPE 
DO 150 II-1.NQD2 
150 INDX( II )«( II-l ) *NNPE+J 

SH( 1 ; NQD2 )-Q8VGATHR( BB( 1 , IT ; NQD2HPE ) , INDX( 1 ; NQD2 ) ; SH( 1 ; NQD2 ) 

I1-(J— 1)*NQD2+1 

B(I1 ,IT;NQD2)-SH(1 ;NQD2) 

15b CONTINUE 

BB( 1 , IT ; NQD21JPE ) - B( 1 , IT; NQD2NPE ) *WTDETEX( 1 ; NQD2NPE ) 

155 CONTINUE 

DO 205 IS-l.HSTR 
DO 205 II-l.NNPE 
NQ-(II-1)*NQD2+1 

205 SIG( NQ , IS ; NQD2) -STRS ( 1 , IS ; NQD2 ) 

DO 400 II-l, NDF 
SUM(1 ;NQD2NPE)-0. 

DO 210 JJ-l.NSTR 
• IT— IBSP (JJ, II) 

213 IF (IT.EQ.O) GO TO 210 

ASSIGN BDESC , BB( 1 , IT ; NQD2NPE ) 

■ SUM( 1 ;NQD2NPE)»SUM( 1 ;NQD2NPE)+BDESC*SIG( 1 , JJ ;NQD2NPE) 

210 CONTINUE 

DO 220 J-l.NNPE 
IJ-(J-1)*3+II 
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NJ-(J-1)*HQD2+1 

220 F0RC(IJ)-Q8SSUM(SUM(NJ ; NQD2) ) 

400 CONTINUE 
RETURN 
END 

. SUBROUTINE MATINV(A,NMAX,N,B, MAX, M, DETERM) 
DIMENSION A( UMAX, UMAX) ,B(NMAX,tlAX) 

DIMENSION IPIVOT(IOO) , INDEX( 100 , 2) ,PIV0T(100) 

IF M-0 IT CALCULATES THE INVERSE ONLY. 

IF M-l IT CALCULATES THE SOL TO AX-B IN B 
INITIALIZATION 

10 DETERM-1.0 
15 DO 20 J-1 ,N 
20 IPIV0T( J)-0 
30 DO 550 1-1 ,N 

SEARCH FOR THE PIVOT ELEMENT 


40 AMAX-0.0 

45 DO 105 J-1,N 

50 IF( IPIVOT( J)-l)60 ,105,60 

60 DO 100 K«1,N 

70 IF(IPIV0T(K)-1)80, 100,740 f 

80 IF( ABS(AMAX)- AT'-' A( J ,K) ) )85 ,100, 100 ! 

85 IRON- J 

90 ICOLUM-K 
95 AMAX- A(J,K) 

100 CONTINUE 

105 CONTINUE 

110 IPIV0T( ICOLUM) "IPIV0T( ICOLUM)+l 

INTERCHANGE ROWS TO PUT ELEIilENG ON DIAGONAL 

i; 

130 IF( IROW-ICOLUH) 140,260,140 ! 

140 DETERM- -DETERM | 

150 DO 200 L-l , N I. 

160 SWAP- A( IROW ,L) J. 

170 A( IROW , L )- A( ICOLUM , L ) j 

200 A(ICOLUM,L)— SWAP } 

205 IF(M)260,260,210 i 

210 DO 250 L-l.M 

220 SWAP- B( IROW ,L) j 

230 B(IROW,L)- B( ICOLUM, L) j 

250 B( ICOLUM, L)- SWAP : 

260 INDEX(I.l)- IRON j 

270 INDEX(I, 2)- ICOLUM l 

310 PIVOT(I)- A< ICOLUM, ICOLUM) 

320 DETERM- DETERM*PIVOT(I) 

DIVIDE PIVOT BY PIVOT ELEMENT 

330 A( ICOLUM, ICOLUlO-l.O 

340 DO 350 L-l ,N 

350 ACICOLUM.L)- A( ICOLUM.L) /PIVOT(I) 

355 IF(M) 380,380,360 

360 DO 370 L-l ,M 

370 B(IC0LUM,L)- B(ICOLUM,L)/PIVOT(I) 

C 

C REDUCE NON-PIVOT ROWS , 

C 

380 DO 550L1-1 ,N 

390 IF ( L1-ICOLUM)400 , 550 , 400 . ' j 

400 T- A( LI, ICOLUM) I 

420 A( LI , ICOLUM)— 0 • 0 f 
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430 DO 450 L-1,N 

450 A(L1,L)- A(L1,L)-A(IC0LUM,L)*T 

455 IF(M) 550,550,460 
460. DO 500 L-1,M 

500 B(Ll.L)- B(L1,L)-B(IC0LUM,L)*T 
550 Continue 
c 

C INTERCHANGE COLUMNS 
C 

600 DO 710 I»1,N 
610 L-tH-l-I 

620 IF(INDEX(L,1)-INDEX(L, 2) >630,710,630 
630 JROH- INDEX(L,1) 

640 JCOLUM- INDEX(L, 2) 

650 DO 705 K«1,N 
660 SWAP- A(K, JROW) 

670 A(K,JROW)- A(K, JCOLUM) 

700 A(K, JCOLUM)- SWAP 
705 CONTINUE 
710 CONTINUE 
740 RETURN 
END 

BLOCK DATA 

COMMON/ GAUSS/ C0RD(8,8) ,WEIGHT(8,8) 

C0MM0N/GENRL/GCR(8 ,3) 

DATA CORD/ 8*0.0, 

A -0.577350269189626,0.577350269189626,6*0.0, 

B -0.774596669241483,0.0,0.774596669241483,5*0.0, 

C -0.861136311594053, -.339981043584856, 0.339981043584856, 

1 0.861136311594053,4*0.0, 

D-0. 906179845938664, -0.538469310105683, 0.0, 0.538469310105683, 

1 0.906179845938664,3*0.0, 

E -0.932469514203152,-0.661209386466265,-0.238619186083197, 

1 +0.238619186083197,0.661209386466265,0.932469514203152,2*0.0, 

F -0.949107912342759, -0.741531185599394,-0.405845151377397,0.0, 
1 0.405845151377397, 0.741531185599394,0.949107912342759,0.0, 

G -0.960289856497536,-0.796666477413627,-0.525532409916329, 

1 -0 . 183434642495650 , 0 . 183434642495650 , 0 . 525532409916329 , 

2 0.796666477413627,0.960289856497536/ 

DATA HEIGHT / 8*0.0, 

A 1.0, 1.0, 6*0.0, 

B 0.555555555555556,0.888888888888889,0.555555555555556,5*0.0, 

C 0.347854845137454,0.652145154862546,0.652145154862546, 

1 0.347854845137454,4*0.0, 

D 0.236926885056189, 0.478628670499366,0.568888888888889, 

1 0.478628670499366, 0.236926885056189,3*0.0, 

E 0.171324492379170,0.360761573048139,0.467913934572691, 

1 0.467913934572691,0.360761573048139,0.171324492379170,2*0.0 , 

F 0.129484966168870,0.279705391489277,0.381830050505119, 

1 0.417959183673469,0.381830050505119,0.279705391489277, 

2 0.129484966168870 ,0.0, 

G 0.101228536290376,0.222381034453374,0.313706645877887, 

1 0.362683783378362, 0.362683783378362, 0.313706645877887, 

2 0.222381034453374,0.101228536290376/ 

DATA GCR/ 0.0,0.0,1.0,1.C,0.0,O.U,1.U,1.0, 

1 1 . 0 , 0 . 0 , 0 . 0 , 1 . 0 , 1 . 0 , 0 . 0 , 0 . 0 , 1 . 0 , 

2 4*0.0, 4*1.0/ 

END 

SUBROUTINE VON(STR,SBAR,PFS) 

C 

*** COMPUTE FLOW VECTOR 

DIMENSION PFS(6),STR(6) 

S1-2*SBAR 

PFS(l)-(2*STR(l)-STR(2)-STRt3))/Sl 
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PFS(2)-(2*STR(2)-STR(1)-STR(3))/S1 

PFS(3)-(2*STR(3)-STR(1)-STR(2))/S1 

PFS(4)-(6*STR(4))/S1 

PFg(5)-(6*STR(5))/Sl 

PFS(6)-(6*STR(6))/S1 

RETURN 

END 

SUBROUTINE DEPL( STR, SBAR, D , HP , DPI) 

DIMENSION STR(&) , STA(6) ,D(6,6) ,DD(6) ,DPL(6,6) 

*** DEP - ELASTIC-PLASTIC MATRIX 

CALL VON ( STR , SBAR, STA) 

BI-D(1,1) 

B2»D(4,4) 

B3-D(l,2) 

DD( 1 )»B1*STA( 1 )+B3* ( STA( 2 )+STA( 3 ) ) 
DD(2)-B1*STA(2)+B3*(STA(1)+STA(3)) 
DD(3)-B1*STA(3)+B3*(STA(1)+STA(2)) 

DD(4)-B2*STA(4) 

DD(5)-B2*STA(5) 

DD(6)«B2*STA(6) 

SD-B1*( STA( 1 ) **2+STA( 2) **2+STA( 3 ) **2 )+2*B3* ( STA( 1 ) *STA( 2 )+ 
C STA( 2 ) *STA( 3 )+STA( 3) * STA( 1 ) )+B2* ( STA( 4)**2+STA( 5) **2+ 

C STA( 6 ) **2 ) 

SD-1.0/(SD+HP) 

DO 10 1-1,6 
DO 10 J-1,6 

10 DPL(I, J)— D(I, J)—DD(I)*DD( J)**SD 

RETURN 
END 

SUBROUTINE SEQU(XT.XXZ) 

C VON MISES YIELD CRITERION. 

DIMENSION XT(6) 

51- O. 5* (XT(1)-XT(2) )**2 

52- 0.5*(XT(2)-XT(3))**2 

53- 0.5*(XT(3)-XT(1))**2 

54- 3*(XT(4)**2) 

55- 3*(XT(5)**2) 

56- 3*(XT(6)**2) 

ST-S1+S2+S3+S4+S5+S6 

XXZ-SQRT(ST) 

RETURN 

END 

SUBROUTINE MULTYS(A,B,N,M,C) 

DIMENSION A(N,M),B(M),C(N) 

DO 10 1-1 ,N 
C(I)-0.0 
DO 10 J-1,11 

10 C(I)-C(I)+A(I,J)*B(J) 

RETURN 

END 

SUBROUTINE PLAS 

COMMON/MAIN/ AA(2000000) ,BB(9600,1) ,D(6,6) ,DINV(6,6) , 
1DISP(80) ,EPS(12400) ,EFEST(1550) .FORCE(IO) ,LINE(100) , 
2L0CAT(10) .LBFOR(IO) ,MB(9600) ,MSUM(9600) ,MPTAB(9600) , 
3MPLAS(12400) ,M0DE(8,1550) ,MPLC(1550) ,N0DX0(700) , 
4N0DY0(700) ,NOD20(700) ,N0DXC(700) ,NODYC(700) ,N0DZC(700) , 
5N0DFIX(80) ,NODLOD(80) ,NDISP(80) ,R(9600) , 

•6SICBAR( 12400) ,SK(24,24) ,T1(9600) ,T2(9600) ,T4(2000) , 
7T3(9600) ,U(3200) ,U0LD(3200) ,V(3200) ,V0LD(3200) ,V2(3200) , 
8W(3200) ,W0LD(3200) ,X(74400) ,XR(3200) ,Y(74400) , 

’9YR(3200) ,Z(9600) ,ZR(3200) 

COMMON/ CNST/ EPS I , SK2 , LMAX , KMAX , DAX , 

1PYLD, SCRIT .YOUNG , POIS , CRACK, PT , WIDTH, P11AX , HP , 

2S BAR, LPRIT , NGAUS , N LAYER , NNODE , , 
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3IN0DX0 , INOD i 0 , INODZO , IHODXC , INODYC , INODZC .LNSTIF , MXNOD , 

4MXNEL .MXGAUS , ICUT , LTOTB , ITNODX.KLU ,NTYP , NLM , 

. 5NDOF, KNEW , N£ P , ERIT , NELM , AM , ROM 
C 

COMMON/MLTNMAT/ YSTRS(20) ,YSTRN(20) ,PLM0DR(20) .NSEGMT 
- COMMON/ D382/NNPE,NDF,NQD,NSTR,NQD2,NNPE2,NQD2NPE,NQD2SR,MXQ2S 
COMMON/VECT/ STRV(8,6) ,STRSV(8,6) ,BMT(64,3) ,WDUM(8) ,XE(8,3) , 

C NCUBE(8) ,DIS(8,3) 

DIMENSION DPL(6 ,6) , AMAT(8 ,3) ,STGAS(6) 

DIMENSION QP(9600) ,STGASV(8,6) ,FFTR(40) 

DIMENSION STR(6) ,02(3200) , INDX(3200) 

DIMENSION YXF.(6) ,ST1(6) ,DPSTRN(6) ,STREP(6) 

DIMENSION PLV(24) 

DATA HALT , GROW/ 4HHALT , 4HGR0W/ 

1000 FORMAT ( 5 X , ' STEP# ' , 12 , 2X , ' TIME : CPU-' , FI 2 . 4 , 2X , 

1 'WALL-' ,F12.3) 

C 

C *** INCREMENT DISPLACEMENTS, FORCES, STRESS & STRAINS TO 1ST YIELD LOAD 
C 

U0LD( 1 ; NNODE ) -U ( 1 ; HNODE ) *PYLD 
VOLD ( 1 ; NNODE ) -V ( 1 ; NNODE ) *P YLD 
WOLD( 1 ; NNODE ) -W ( 1 ; NNODE )*PYLD 
R(1 ;NDOF)-R( 1 ; NDOF)*PYLD 
X( 1 ; MXQ2S )-X( 1 ; MXQ2S )*PYLD 
Y(1 ;MXQ2S)-Y(1 ;MXQ2S)*PYLD 
C *** ZEROING 

QP(l;ND0F)-0.0 
EPS( l;MXGAUS)-0 .0 
MPLAS ( 1 ; MXGAU S ) — 0 
C 

C *** READ DATA 

C 

READ(5,11) PCT,ERIT,MAXIT ,N0DE1 ,N0DE2 ,NELE1 ,NELE2 PLAS 

11 FORMAT (2E1 0.3, 515) PLAS 

READ( 5,312) P.WORD 

312 FORMAT(E10.3,1X,A4) 

WRITE( 6 , 313 ) P , PCT , ERIT , MAXIT , N0DE1 , N0DE2 , NELE1 , NELE2 

313 FORMAT(///10X, 'TOTAL LOAD FACTOR-' ,F10. 4 

l /10X,' INCREMENTAL LOAD FACTOR-' ,F10. 4/ 10X,' ALLOWABLE ERROR ON STRE 
2SS-',F10.4/10X, 'MAXIMUM NUMBER OF ITERATION-' ,14 
3/10X, 'PRINT DISPLACEMENTS AT NODES', 15,' TO', 15 
4/10X, 'PRINT STRESSES IN ELEMENTS' ,15,' TO', 15) 

CALL Q3CL0CKS( CPU, WALL) 

IZIP1-0 

WRITE(6 , 1000) IZIP1, CPU, WALL 
PT-PYLD 

CALL PLOUT(NODE1 ,N0DE2 ,NELE1 ,NELE2) 

CALL Q3CL0CKSC CPU, WALL) 

IZIP1-1 

WRITE(6 , 1000) IZIP1, CPU, WALL 
IF(NEP.EQ.O) STOP 
DELP-PCT*PYLD 
NPL-0 

20 PMAX-PT 

PT-PT+DELP 

NPLOT-O 

IF(PT.GE.P . AND.DELP.GT. 0.0) GO TO 25 
IF(PT.LE.P. AND.DELP.LT. 0.0) GO TO 25 
GO TO 26 

25 - PT-P 

NPLQT-1 

26 ’ DELPO-PT-PMAX 

NPL-NPL+1 

KLU-0 

ICON-O 

V2( 1 ;NNODE)-VOLD( 1; NNODE) 


c 

C HOLDING APPLIED LOAD CONSTANT -ITERATE UNTIL SOLUTION CONVERGE 

C 

PTOY-PT/PYLD 
NBREAK-0 
• GO TO '45 

35 NL— 1 

36 NL-NL+1 
NPLOT-O 

IF(NL.GT.NLM) GO TO 91 
ANL-NL 

DO 133 JIS-l.ICUT 

FFTR( JIS)-FORCE( JIS)*( 1 .-ANL/NLM) 

133 URITE(6 , 167)FFTR( JIS) ,HL 

167 F0RMAT(10X, 'CRACK-TIP FORCE-' ,E16 . 7 , 'AT STEP', 12) 

45 DO 50 ITER-l.MAXIT 

tlC-0 

65 BB( 1,1; NDOF)-QP( 1 ; ND0F)+R( 1 ;NDOF) *PTOY 

IF(KLU.EQ.l) GO TO 830 
GO TO 831 

830 DO 134 JIT-1, ICUT 
NOM-LOCAT(JIT) 

NFL«3*N0M-1 

134 BB ( NFL , 1 )- BB ( NFL , 1 )+FFTR( JIT ) 

831 CONTINUE 
IF AC-1 

NC-1 

CALL S Y11BAN(LNSTIF , NDOF , MB ,MSUM , AA, 1 , BB , IFAC , T1 , IERR , 

1 ALP , Z , T2 , T3 , T4 , NC) 

DO 70 N-l.NNODE 
70 INDX(N)-(N-1)*3+1 

U2( 1 ;NNODE)- Q8VGATHRC BB( 1 , 1;ND0F) ,INDX( 1 ; NNODE) ; U2( 1 ; NNODE ) ) 

U ( 1 ; NNODE ) -U2 ( 1 ; NNODE ) -UOLD ( 1 ; NNODE ) 

UOLD( 1 ;NN0DE)-U2( 1 ;NN0DE) 

INDX ( 1 ; NNODE) -INDX( 1 ; NN0DE)+1 

U2(l; NNODE)- Q8VGATHR( BB( 1,1; NDOF) , INDX( 1 ; NNODE ) ; U2 ( 1 ; NNODE ) ) 

V ( 1 ; NNODE ) -U2 ( 1 ; NNODE )-VOLD ( 1 ; NNODE ) 

VOLD( 1 ; NN0DE)-U2 ( 1 ; NNODE ) 

INDX( 1 ; NNODE ) -INDXQ ; NN0DE)+1 

U2( 1 ;NNODE)- Q8VGATHR( BB( 1,1; NDOF) , INDX( 1 ; NNODE) ; U2( 1 ; NNODE) ) 

W ( 1 ; NNODE ) -U2 ( 1 ; NNODE )-WOLD ( 1 ; NNODE ) 

WOLD( 1 ; NNODE )-U2( 1 ; NNODE) 

C 

C COMPUTE TOTAL STRAIN INCREMENTS FROM DISPLACEMENT INCEMENTS. 

C COMPUTE ELASTIC STRESS INCREMENTS AND ADD TO CURRENT STRESSES. 

C CHECK YIELD CONDITION FOR PLASTIC ELEMENTS. 

C 

IGAUSP-0 
DO 80 I-l.NELM 
DO 75 J-1,8 
NCUBE( J)-MODE( J ,1) 

Nl-NCUBE(J) 

DIS( J, l)-U(Nl) 

DIS( J , 2)-V (Nl) 

75 DISC J,3)“W(N1) 

C *** 

CALL CORDIN(NCUBE,MXNOD,XR,YR, ZR.XE) 

CALL -STRESS ( DIS , XE , D , STRV , STRSV , BMT , WDUM) 

STGASV( 1 , 1 ;NQD2SR)-0 .0 
• IL0C-(I-1)*NQD2SR 
DO 76 IG-1 , NQD2 
IGAUSP-IGAUSP+1 
' SBAR-SIGBAR(IGAUSP) 

DO 77 JS-1,6 

STR( JS)-STRSV(IG,JS) 

YTR(JS)-STRV(IG, JS) 
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JS1«IL0C+IG+NQD2*(JS-1) 

77 ST1( JS)-X( JS1) 

CALL SEQU(ST1,S1) 

AU-STR(l) 

A22-STR(2) 

' A33-STR(3) 

A44-STR(4) 

A55»STR(5) 

A66-STR(6) 

su-sti(i) 

S22-ST1(2) 

S33-ST1C3) 

S44-ST1(4) 

S55-ST1(5) 

S66-ST1(6) 

ST1(1; 6)-STl(l; 6)+STR(l; 6) 

CALL SEQU(ST1,S2) 

C *** IF SOL CONVERGED THEN GOTO 801 
IF(ICON.EQ.l) GO TO 801 
IF(S2.LT.S1) MPLAS(IGAUSP)-0 
IF(S2.LT.S1) GO TO 801 

IF(M?LAS(IGAUSP) .WE.O.AND.ITER.GT.l) GO TO 74 
IF(S2.LE.SBAR) GO TO 801 
MPLAS ( IGAU SP ) -IGAUSP 
74 MP-1 

C CHECK FOR CONVERGENCE. 

IF(ABS(S2-SBAR) .GT.ERIT) MC-1 

A-Al l**2+A22**2+A33**2+3* ( A44**2 )+3* ( A55**2 )+3* ( A66**2 ) 

1 -(A11*A22)-(A22*A33)-(A33*A11) 

B22»S11*(2*A11-A22-A33)+S22*(2*A22-All-A33)+S33*(2*A33-A22-All) 
1 +6*S44*A44+6*S55*A55+6*S66*A66 
C»S1**2-SBAR**2 

IF(A.LT.EPSI) GO TO 8 
IF(ITER.EQ.2) GO TO 200 
DELTA“B22**2-4*A*C 
IF (DELTA) 2 00 ,40, 40 
200 PX-(SBAR-S1)/(S2-S1) 

GO TO 231 

40 P0NE-(-B22+SQRT(DELTA))/(2.*A) 

PTWO- (-B22-SQRT( DELTA) ) / ( 2.* A) 

PX-PONE 

IF ( ABS ( PC NE ) . GT . AB S ( PTOO ) ) PX-PTWO 
231 CONTINUE 

PXD-l.-PX 

YTR(1; 6)-YTR(l; 6)*PXD 
STR( 1 ; 6) -STR( 1 ; 6) *PXD 

IF(ROM.LE.O. .AND. NSEQ1T.EQ.0) HP-A11*Y0UNG 

IF(ROM.LE.O. .AND. NSEG11T.GT.0) HP“FNMA1’(SBAR, YOUNG, EPS(IGAUSP) , 
CIGAUSP) 

IF( ROM . GT . 0 . ) HP-ROM**AM*SBAR** ( 1 . -AM) /AM 
CALL DEPL(ST1 , SBAR,D,HP,DPL) 

CALL MULTYS(DPL,YTR,6,6,STREP) 

CALL MULTYS(DINV, STREP, 6, 6 .DPSTRN) 

DPSTRN(1 ; 6)-YTR(l ; 6)-DPSTRN( 1 ; 6) 

STGAS (1;6)-STR(1; 6 )-STREP ( 1 ; 6 ) 

CALL ERTA( DPSTRN, SMA) 

EPS(IGAUSP)-EPS(IGAUSP)+SMA 

HC-1.0 

SIGBAR( IGAUSP)-SBAR+HC*HP*SMA 
■ DO 12 IP-1,6 
12 STGASV(IG,IP)-STGAS(IP) 

8 CONTINUE 

C ' IF(I.EQ.l) WRITE(6 ,400) SBAR,EPS( IGAUSP) , SIGBAR( IGAUSP) , HP, SI ,S2 
C400 F0RMAT(5X, 'CONSTANTS: '.6F12.3) 

7b CONTINUE 

STRSV( 1 , 1 ;NQD2SR)-STRSV(1 , 1;NQD2SR)-STGASV(1 , 1 ;NQD2SR) 



CALL FORCEP ( BMT , HDUM , STGAS V , PL V) 

DO 455 XT-1,8 
XX-3* IT- 2 
Nl-NCUBE(IT) 

NU«3*Nl-2 

. QP(NU) -QP(NU)+PLV( IX) 

QP(NU+1) -QP ( NU+1 )+PL V ( IX+1 ) ' 

455 QP(NU+2)-QP(NU+2)+PLV(IX+2) 

801 IL0C1-IL0C+1 

X( IL0C1 ; NQD25R)-X( IL0C1 ; NQD2SR)+STRSV (1,1; NQD2SR) 

Y( IL0C1 ; NQD2SR)-Y( IL0C1 ; NQD2SR)+STRV( 1,1; NQD2SR) 

80 CONTINUE 

IF( ICON. EQ. 1) GO TO 90 
IF(MP.EQ.O) GO TO 90 
IF(MCcEQ. 1) GO TO 49 
WRITE (6, 52) ITER 

52 F0RMAT(10X, 'SOLUTION CONVERGED IN' ,14 .'ITERATIONS' ) 

ICON-1 

GO TO 49 

53 WRITE (6,54) ITER 

54 FORMAT(10X,'N0 CONVERGENCE IN ' ,14, 'ITERATION') 

CALL PLOUT( N0DE1 , N0DE2 , HELE1 , NELE2) 

GO TO 999 

49 IF(ITER.EQ.MAXIT) GO TO 53 
IF(KLU.EQ. 1 .AND.NL.EQ.O) GO TO 90 

50 CONTINUE 

90 CONTINUE 
MP-0 
ICON-O 

IF(KLU.EQ.l) GO TO 36 

91 CONTINUE 
CALL CONTACT 
IF(KNEW.EQ.l) GO TO 45 

IF(NPLOT.EQ.l) CALL PLOUT(NODEl ,N0DE2 .NELEl ,NELE2) 
IF(NPL.EQ.NEP) CALL PLOUT ( NODE 1 , N0DE2 , NELE 1 , NELE2 ) 
IF(NPL.EQ.NEP) NPL-0 
IF(NTYP.EQ.l) CALL BREAK 
IF(NTYP.EQ.l) GO TO 100 • 

IF(KLU.EQ.l) GO TO 1C, 2 

IF(NPLOT.EQ. 1. AND. WORD. EQ. GROW) CALL BREAK 
100 CONTINUE 

IF(KLU.EQ.2) GO TO 999 
IF(KLU.EQ.3) GO TO 999 
IF(KLU.EQ.l) GO TO 35 
IF(NPLQT.EQ.l) GO TO 99 
GO TO 20 
99 CONTINUE 

102 DELP-- DELP 

READ(5,121)P,W0RD 
121 FORMAT(E10.3,1X,A4) 

IF(WORD.EQ.HALT) GO TO 999 
NPL-0 

MPLAS( 1 ; MXGAUS )-0 
GO TO 20 
999 RETURN 

END 

SUBROUTINE 3RTA( A , B) 

DIMENSION A(6) 

X-A(l) 

Y«A(2) 

Z-A(3) 

XY-A(4)/2 . 

YZ-A(5)/2. 

ZX-A(6)/2. 

51- (X-Y)**2 

52- (Y-Z)**2 
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S3«(Z-X)**2 

54- 6*(XY)**2 

55- 6*(YZ)**2 

56- 6*(ZX)**2 

STOT- S 1+S2+S3+S4+S5+S6 
' SPO-SQRT(STOT) 

SO-SQRT(2.)/3. 

B-SO*SPO 

RETURN 

END 

SUBROUTINE PLOUT( NODE1 , NODE2 , NELE1 , NELE2 ) 

COMMON/ MAIN/ AA(2000000) ,BB(9600,1) ,D(6,6) ,DINV(6,6) , 

1DISP(80) ,EPS(12400) ,EFEST(L550) ,F0RCE(10) ,LINE(100) , 

2LOCAT(10) ,LBFOR(10) ,MB(9600) ,MSUM(9600) ,MPTAB(9600) , 

3MPLAS( 12400) ,M0DE(8 , 1550) ,MPLC(1550) ,NODXO(700) , 

4N0DY0(700) ,NODZO(700) ,N0DXC(700) ,N0DYC(700) ,NODZC(700) , 
5N0DFIX(80) ,N0DL0D(80) ,NDISP(80) ,R(9600) , 

6SIGBAR( 12400) ,SK(24,24) ,T1(9600) ,T2(9600) ,T4(2000) , 

7T3(9600) ,U(3200) ,UOLD(3200) ,V(3200) ,V0LD(3200) ,V2(3200) , 

8W( 3200) ,WOLD( 3200) ,X(74400) ,XR(3200) ,Y(74400) , 

9YR(3200) ,Z(9600) ,ZR(3200) 

COMMON/ CNST/EPS I , SK2 , U-1AX , KMAX ,DAX , 

1PYLD , SCRIT , YOUNG , POIS , CRACK, PT , WIDTH , PMAX , HP , 

2SBAR , LPRIT , NGAUS , NLAYER , NNODE , 

31NODXO , INODYO , 1N0DZ0 , INODXC , INODYC , INODZC.LNSTIF , MXNOD , 
4MXNEL.MXGAUS ,ICUT , LTOTB , ITNODX, KLU ,NTYP ,NIM, 

5ND0F, KNEW , NEP , ER1T , NELM , AM .ROM 

C0MM0N/D382/NNPE.NDF ,NQD,NSTR,NQD2,NNPE2 ,NQD2NPE,NQD2SR,MXQ2S 
COMMON/VECT/ STRV(8,6) ,STRSV(8,6) ,BMT(64,3) ,WDUM(8) ,XE(8,3), 

C NCUBE(8) ,DIS(8,3) 

DIMENSION STRS(6) ,FRC(24) ,FORCEX(3200) ,FORCEY(3200) ,F0RCEZ(3200) 
*'** OUTPUT ROUTINE 

WRITE(6 , 10)PT, CRACK, WIDTH 

10 FORMATC/.10X, 'APPLIED LOAD-' ,E12. 5, 8X, 'CRACK-' , 

1 F10. 5 , IOX, 'WIDTH-' , F10, 5/) 

IF( CRACK. LT.EPSI) GO TO 20 
WRITE(6,15) 

1 5 FORMAT ( 12X, 'NODE' , 5X, 'X' , IOX, ' Y' , IOX , ' V , IOX, 

1 'COD'/) 

CRACK1-CRACK+EPSI 
CRACK2«CRACK-10*DAX 
DO 16 1-1, INODYO 
L-NODYO(I) 

IF(XR(L) .LT.CRACK2) GO TO 16 
IF(YR(L) .GT.EPSI) GO TO 16 
IF(XR(L) .GT.CRACK1) GO TO 16 
WRITE(6,25)L,XR(L) ,YR(L) ,ZR(L) ,VOLD(L) 

25 F0RMATC5X, 14 , 4E14 . 6) 

16 CONTINUE 

20 CONTINUE 

WRITE(6,30) 

30 FORMAT( // , 30X, 'DISPLACEMENTS'/ /12X,' NODE' , 14X,'U' , 

1 18X,'V',18X,'W') 

DO 12 N-N0DE1 ,N0DE2 

12 WRITE(6 ,22) «,UOLD(N) ,VOLD(N) ,UOLD(N) 

22 FORMAT (IOX, 15, 5X,3(E13.6, 3X) ) 

WRITE(6,35) 

35 F0RMAT(//20X,'STRESSES AND STRAINS' , IOX, 1H* ,3X, 

1 'DENOTES PLASTIC ELEMENTS' ,30X, 'EFFECTIVE'/ / , 5X> 

' 1 'ELEMENT' , 5X , ' S IGX' , 8X , ' S IGY' , 8X, ' S IGZ ' , 8X , ' TAUXY' , 

1 8X, 'TAUYZ' ,8X, 'TAUZX' ,8X, 'EFFE-STRESS'/ /) 

C 

C *** CALCULATE NODAL FORCES AND GAUSS POINT STRESSES 
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WRITE(6 , 25) tIGAUS 
NGAUS-2 

| FORCEX(1;NNODE)-0.0 
FORCEY ( 1 ; NNODE) -0 . 0 
- F ORCEZ ( 1 ; NNODE )«0.0 
N-0 

IGAUSP-0 

DO 300 IE-1, HELM 
NGUBE(1;8)-M0DE(1,IE;8) 

CALL C0RDIN(NCUBE,MXN0D,XR,YK 
IL0C1- ( IE-1 ) *NQD2SR+1 
SXRSV( 1 , 1; NQD2SR)- X(IL0C1;NC 
CALL CDER( XE , NGAUS , BMT , WDUM , 2 
CALL FORCEP(BHT,WDUtl, STRSV.FR 
DO 352 1-1,8 

Il-NCUBE(I) 

IX-l*3-2 

FORCEX( I 1 )-FORCEX( I1)+FRC( IX) 
FORCE Y ( I 1 ) -FORCE Y ( II )+FRC ( IX+ 
352 FORCEZ( Il)-PCkCEZ( I1)+FRC( IX+ 
DO 350 IG-1.NQD2 
IGAUSP-IGAUSP+1 
DO 351 1-1,6 

351 STRS ( I ) - STRSV ( IG , I ) 

SIGO- ( 1-ERIT ) *SIGBAR( IGAUSP ) 
CALL SEQU ( STRS , STP ) 
IF(STP.GE.SIGO) GO TO 42 
WRITE(6 ,43)IE, IG. ( STRS(L) ,L= 
43 FORMAT (5X, 16, I2.5X, 7E12. 5) 

GO TO 350 

42 WRITE(6 , 45) IE, IG, ( STRS(L) ,L= 

45 F0RMAT(4X, 1H* , 16 , I2,5X,7E12 

XGAUSP-KGAUSP+1 
350 CONTINUE 

IF(KGAUSP.LE.O) GOTO 300 
93 N-N+l 

HPLC(N)— IE 
300 CONTINUE 

NOPL-N 
C 

C *** 

C 

WRITE(6 ,371) 

DO 370 IN-N0DE1 ,N0DE2 

370 URITE(6 ,372) IN.FORCEX(IN) ,FC 

371 FORMAT(1H1///10X,'NODE #',5X, 

372 FORMAT( 10X,I5,2X,3(E12.4,2X) ) 
C 

C *** 

C 

WRITE( 6 , 380) 

380 F0RMAT(//10X,'L1ST OF PLASTI 
URITE(6 ,381) (MPLC(I) ,1-1, N) 

381 F0RMAT(5X, 2015) 

997 RETURN 

END 

SUBROUTINE CONTACT 
COMMON/MAIN/ AA(2000000) ,BB(96 
•1DISP(80) ,EPS(12400) ,EFEST(155 
2L0CAT( 10) ,LBF0R(10) ,MB(9600) , 
' 3MPLAS( 12400) ,M0DE(8 , 1550) ,MP1 
4N0DYO(70O) ,NODZ0(700) ,N0DXC(7 
5N0DFIX(80) ,NODLOD(80) ,NDISP(8 
6SIGBAR( 12400) ,SK(24 ,24) ,T1(9( 
7T3(9600) ,0(3200) ,UOLD(3200) ,\ 


l.XE) 

3R) 


6) ,STP 


6) ,STP 

) 


EY(IN) ,F0RCEZ(IN) 

ORCEX' , 8X , 'FORCEY' , 7X, ' FORCEZ' ) 


ELEMENTS'//) 


,1) ,D(6,6) ,DINV(6,6) , 
.FORCE (10) ,LINE( 100) , 
■UM(9600) ,MPTAB(9600) , 
1550) ,N0DXO(700) , 

) ,N0DYC(700) ,NODZC(700) , 
,R(9600), 

0 ,T2(9600) ,14(2000) , 

200) , V0LD(3200) ,V2(3200) , 



SWC3200) ,W0LD(3200) ,X(74400) ,XR(3200) ;Y(74400) 

9YR(3200) ,Z(9600) ,ZR(3200) 

. COMMON/ CNST/ EPSI, SK2 ,LMAX,KMAX,DAX, 

. 1PYLD, SCRIT, YOUNG, PO IS .CRACK, FT, WIDTH, PMAX.'tlP , 

2SBAR , LPRIT , NGAUS , NLAYER , NNODE , 

3 IflODXO , INODYO , INODZO , INODXC , INODYC , INODZC .LNSTIF , MXNOD , 
4MXNEL , MXGAUS , ICUT , LTOTB , ITNODX, KLU , NTYP , NLM , 

5NDOF , KNEW , HEP , ERIT , NELM, AM, ROM 

C CHANGES SPRING STIFNESS IF CRACK CLOSES OR OPENS. 

C DONT FORGET TO PUT THE COMMON HERE. 

WRITE (6, 95) 

95 FORMAT(5X, 'CALLING FROM THE CONTACT') 

DO 10 I- 1, IliODYO 
L-NODYO(I) 

KNEW-0 

CX«XR(L)-CRACK 
IF(CX.LT.-EPSI) GO TO 9 
GO TO 40 

9 MV-3*L-1 
MPTX-MPTAB(HV) 

IF(VOLD(L) .LE.O.O) MPTAB(NV)-1 
IF( VOLD(L) .GT.O. 0) MPTAB(NV)-0 
IF(HPTX.NE.MPTAB(NV) ) KNEW-1 
IF(KNEW.EQ.O) GO TO 40 

PN-PT-( ( PT-PMAX) / ( VOLD(L )-V2 ( L) ) ) * VOLD(L) 

Z(NV)“1 . 0 

NC-NV 

ALP-SK2 

IF(MPTAB(NV) .EQ.O) ALP— SK2 
IFAC-3 

CALL S YMBAN(LNSTIF , NDOF , MB , MSUM.AA, 1 , BB , IFAC , T1 , IERR, 

1 ALP,Z ,T2 ,T3 ,T4 ,NC) 

IF(MPTAB(NV) .EQ.O) . WRITE(6 , 20)L,PN 
20 F0RMAT(/,2X, 'NODE', 13, 'OPENED AT',F8.3) 

IF(MPTAB(HV).EQ.l) WRITE<6 ,30)L,PN 
30 F0RI1AT(/,2X, 'NODE', 13, 'CLOSED AT',F8.3) 

40 CONTINUE 

10 CONTINUE 
RETURN 
END 

SUBROUTINE BREAK 

COMMON/MAIN/ AA( 2000000) ,BB(9600,1) ,D(6,6) ,DINV(6,6) , 
1DISP(80) ,EPS(12400) ,EFEST(1550) ,F0RCE(10) ,LINE(100) , 
2L0CAT( 10) .LBFOR(IO) ,ilB(9600) ,HSUM(9600) ,MPTAB(9600) , 
3MPLAS( 12400) ,H0DE(8 , 1550) , MPLC( 1550) ,NODXO(700) , 
4N0DYO(7O0) ,NODZO(700) ,N0DXC(700) ,NODYC(700) ,N0DZC(700) , 
5NODFIX(80) ,N0DL0D(80) ,NDISP(80) ,R(9600) , 

6SIGBAR( 12400) ,SK(24,24) ,T1(9600) ,T2(9600) ,T4(2000) , 
7T3(9600) ,U(3200) ,UOLD(3200) ,V(3200) ,V0LD(3200) ,V2(3200) , 
8W(3200) ,WOLD(3200) ,X(74400) ,XR(3200) ,Y(74400) , 

9YR(3200) ,Z(9600) ,ZR(3200) 

COMMON/ CNST/ EPS I , SK2 , LMAX , KMAX , DAX , 

1PYLD, SCRIT .YOUNG , POIS , CRACK, PT , WIDTH , FMAX , HP , 

2S BAR , L PR IT , NGAUS , NLAYER , NNODE , 

3IN0DX0 , INODYO , INODZO , INODXC , INODYC , INODZC .LNSTIF , MXNOD , 
4MXNEL , MXGAUS , ICUT , LTOTB , ITNODX, KLU , NTYP , NLM , 

5ND0F , KNEW , NE P , ERIT , NELM , AM , ROM 

C COMMON GOES HERE. 

DO 8 I- 1, INODYO 
L-NODYO(I) 

C2“ABS( ZR(L)-ZCOR) 

IF(C2.LT.EPSI) GO TO 10 
GO TO 8 

10 CX«ABS(XR(L)-CRACK) 
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IF(CX.LT.EPSI) GO TO 9 

8 CONTINUE 

9 . IF(NTYP.EQ.O) GO TO 12 

JSUM-O 

DO 90 10-1, INODYO 
NS-'NODYO(IO) 

C1-ABS(ZR(NS)-ZC0R) 

IF(Cl.LT.EPSI) GO TO 130 
GO TO 90 

130 JSUM-JSUMfl 

LINE(JSUM)“NS 

90 CONTINUE 
ITNODX-JSUM 
DIST-YOUNG 

DO 91 IP-l.ITNODX 
LM-LINE(IP) 

CSD“XR(L)-XR(LM) 

IF(CSD.LE.O.O) GO TO 91 

IF( CSD . GT . EPSI . AND. CSD. LT. DIST) DIST-CSD 

91 CONTINUE 
ICAM-0 
MODF-LM 

DO 92 LU“1 , INODYO 
C10«ABS(XR(LU)-XR(M0DF)) 

IF(CIO.LT.EPSI) GO TO 96 
GO TO 92 

96 ICAM-ICA1H-1 

92 LBFOR(ICAM)-LU 
LTOTB-ICAM 

DO 94 JO-l.LTOTB 
LA-LBFOR(JO) 

( STR-VOLD(LA) 

WRITE(6,15)LA,STR,PT 

15 F0RMAT(2X, 'CRACK-TIP NODE' ,14,'HAD' ,E11.4,'CT0D AT' 

1 ,E11.4) 

KLU-0 

IF(STR.GE.(0.98*SCRIT))KLU-1 
IF(KLU.EQ.l) NBREAK-NBREAK+1 
IF(KLU.Eq.O) GO TO 997 
94 CONTINUE 

12 II-O 

DO 13 JJ-1, INODYO 
LL-NODYO(JJ) 

C3-XR(LL) 

Cl-XR(L) 

C5-YR(L) 

C6-YR(LL) 

C7“ABS(C5-C6) 

C4“ABS(C3-C1) 

IF( C4.LT.EPSI.AND.C7.LT. EPSI) GO TO 155 
GO TO 13 
155 II-II+l 

LOCAT(II)-LL 

13 CONTINUE 
ICUT-II 

DO 16 JI-l.ICUT 
LC-LOCAT(JI) 

NV“3*LC-1 

MPTAB(NV)-0 

KLU-1 

IF(NBREAX.EQ.5) KLU-3 

ALP— SK2 

Z(NV)-1. 

NC-NV 

IFAC-3 

CALL SYI1BAN ( LNSTIF , NDOF , MB , MSUM , AA , 1 , BB , IF AC , T1 , IERR , 
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1 ALP.Z.TZ.TS^.NC) 

WRITE(6,100)LC,PT 

100 F0RMAT(/,2X, 'NODE' ,14, 'BROKE AT',F8.3/) 

16 ' CONTINUE 

• CMIN-1.0E+10 

DO 36 ’ JS“1 ,ICUT 
LO-LOCAT(JS) 

36 FORCE(JS)— SK2*VOLD(LO) 

DO 60 IJ-l.INODYO 
LT-NODYO(IJ) 

C5“XR(LT)-XR(L) 

60 1F( C5.GT.EPSI.AND.C5.LT.CMIN) CMIN-C5 

CRACK-CRACK+CUIN 

997 RETURN 
END 



SUBROUTINE SYMBANQ1AXN , N , M.MSUM , A , NRHS , B , IF AC , P , IERR, ALP , Z , W , D , T4 , 

SYMBAN 

2 

c 

C NC) 

SYMBAN 

3 

c 

SOLVE MATRIX EQUATION AX-B WHERE A IS SYMMETRIC POSITIVE DEFINITE 

SYMBAN - 

4 

c 

AND B IS A MATRIX OF CONSTANT VECTORS. 

SYMBAN 

5 

c 

SYMMETRY AND BAND WIDTH OF MATRIX A IS ACCOUNTED FOR BY STORING A(I 

SYMBAN 

6 

c 

IN LOWER TRIANGULAR MATRIX (INCLUDING DIAGONAL)- BY ROWS 

SYMBAN 

7 

c 


SYMBAN 

8 

c 

MAXN - MAXIMUM DIMENSION OF MATRIX A 

SYMBAN 

9 

c 


SYMBAN 

10 

c 

N - NUMBER OF ROWS OF MATRIX A 

SYMBAN 

11 

c 


SYMBAN 

12 

c 

M(T) - BAND WIDTHS (LARGEST NUMBER OF NON-ZERO COLUMNS TO THE LEFT 

SYMBAN 

13 

c 

AND INCLUDING THE DIAGONAL OF ROW I IN MATRIX A 

SYMBAN 

14 

c 

. M(I) MUST BE GREATER THAN OR EQUAL TO 2 WITH M(l)-1 

SYMBAN 

15 

c 


SYMBAN 

16 

c 

MSUM(I) - AN ARRAY COMPUTED IN SYMBAN 

SYMBAN 

17 

c 


SYMBAN 

18 

c 

c 

c 

c 

c 

T4 ( MXBND ) WORKING STORAGE OF LENGH MAX BAND WIDTH. 



MXBND - MAX BAND WIDTH. T4 IS DIMENSIONED ONLY IN MAIN'. 
NRHS - NUMBER OF RIGHT-HAND SIDES OF COLUMN VECTOR B 

SYMBAN 

19 

c 


SYMBAN 

20 

c 

IF AC - INPUT INTEGER SPECIFYING WHETHER OR NOT A CHOLESKY 

SYMBAN 

21 

c 

DECOMPOSITION OF MATRIX A IS TO BE COMPUTED 

SYMBAN 

22 

c 

WHERE A - L*D*L**T 

SYMBAN 

23 

c 


SYMBAN 

24 

c 

- 0 CHOLESKY DECOMPOSITION IS COMPUTED. IFAC SET TO 1. 

SYMBAN 

25 

c 


SYMBAN 

26 

c 

- 1 THE CHOLESKY DECOMPOSED FORM OF MATRIX A IS INPUT. 

SYMBAN 

27 

c 

SOLUTION IS RETURNED IN B. 

SYMBAN 

28 

c 


SYMBAN 

29 

c 

- 2 THE CHOLESKY DECOMPOSITION OF MATRIX A IS MODIFIED. 

SYMBAN 

30 

c 

MODIFICATION IS OF THE FORM A - A + ALP*Z*Z**T. 

SYMBAN 

31 

c 

NC IS THE FIRST NONZERO ROW IN COLUMN VECTOR Z. 

SYMBAN 

32 

c 

IFAC IS SET TO 1 AND Z SET TO 0 UPON RETURN. 

SYMBAN 

33 

c 

SOLUTION IS RETURNED IN B. 

SYMBAN 

34 

c 


SYMBAN 

35 

c 

- 3 ONLY THE CHOLESKY DECOMPOSITION OF MATRIX A IS MODIFIED 

SYMBAN 

36 

c 

IFAC IS SET TO 1 AND Z SET TO 0 UPON RETURN. 

SYMBAN 

37 

c 


SYMBAN 

38 


DIMENSION A(MAXN) ,B(N,NRHS) ,P(N) ,W(N) ,D(N) ,Z(N) ,M(N) ,MSUM(N) ,T4(1) 

SYMBAN 

39 


IF (IFAC.GT.O) GO TO 11 

NEW 

180 


LL-1 

SYMBAN 

41 


DO 10 1-1 ,N 

SYMBAN 

42 


' MSUM(I)"LL 

NEW 

181 


. ' LL“LL+M(I) 

SYMBAN 

44 


10 CONTINUE 

SYMBAN 

45 


11 CONTINUE 

NEW 

182 


IERR-0 
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CALL QSOLV( A , MSUU , M , B , P , N , IFAC , D , W , Z , T4 , ALP , NC , IERR) HEW 183 

RETURN SYMBAN 48 

, END SYMBAN 49 

. SUBROUTINE QSOLV(AR,IB,IL,B,DI,N,NFACT,D,W,Z,T,ALP,NC,IERR) 

- DIMENSION AR(1),IB(1), IL(1),B(1), DI(1), D(l) ,W(1) ,Z(1) NEW 

. DIMENSION T( 1) NEW 

DESCRIPTOR AV.BV QSOLV 

IF(NFACT .GE.2) GO TO 300 NEW 

IF (NFACT.NE.O) GO TO 160 QSOLV 

C FACTOR QSOLV 

DO 100 I»i,N QSOLV 

ICI“I-IL( I)+l QSOLV 

T(1;IL(I))«AR(IB(I);IL(I)) QSOLV 

N1“IB(I)+I— ICI QSOLV 

AR(N1)— 1 QSOLV 

DO 100 J-ICI.I QSOLV 

ICJ“J-IL(J)+1 QSOLV 

KS“MAXO(ICI,ICJ) QSOLV 

Nl-KS-ICI+1 QSOLV 

N2-J-KS+1 QSOLV 

ASSIGN AV,T(N1;N2) QSOLV 

N1-IB(J)+KS-ICJ QSOLV 

ASSIGN BV,AR(N1;N2) QSOLV 

C“Q8S DOT ( AV , BV) QSOLV 

Nl-J-ICI+1. QSOLV 

T(N1)— C QSOLV 

IF (J.EQ.I) GO TO 110 QSOLV 

N2-IB(I)+J— ICI QSOLV 

AR(N2)-T(N1)*DI(J) QSOLV 

GO TO 100 . QSOLV 

110 CONTINUE QSOLV 

IF(T(N1) .LE.O.O) GOTO 999 

DI(I)-1/T(N1) QSOLV 

D(I)-T(N1) 

100 CONTINUE QSOLV 

C FORWARD SUBSTITUTION QSOLV 

160 CONTINUE QSOLV 

DO 200 I”1 ,N QSOLV 

ICI«I-IL(I)+1 QSOLV 

ASSIGN AV, AR(IB(I) ; IL(I) ) QSOLV 

ASSIGN BV,B(ICI;IL(I)) QSOLV 

C-Q8SD0T(AV,BV) QSOLV 

B(I)— C QSOLV 

200 CONTINUE QSOLV 

C DIAGONAL QSOLV 

B(1;N)-B(1;N)*DI(1;N) QSOLV 

C BACKWARD SUBSTITUTION QSOLV 

NM1-N-1 QSOLV 

DO 400 II“1,NM1 QSOLV 

I-N-II+1 QSOLV 

ICI-I-IL(I)+1 QS0LV 

IF (ICI.GE.I) GO TO 400 QSOLV 

B( ICI ; I-ICI)«B( ICI; I-ICI)-AR( IB( I) ; I-ICI)*B( I ) QSOLV 

400 CONTINUE QSOLV 

C SOLUTION IS HOW IN B QSOLV 

C D(1;N)-1.0/DI(1;N) NEW 

RETURN QSOLV 

0 ' NEW 

c********************* HFACT ■ 2 OR 3 ******************* ************ 

C NEW 

300. CONTINUE NEW 

' ' DO 310 J-NC.N NEW 

. ’ D(J)«D(J)+AI.P*Z(J)*Z(J) NEW 

BETA«ALP*Z(J)/D(J) NEW 

ALP-ALP/ (D(J)*DI(J)) NEW 

DI(J)«1.0/D(J) NEW 
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IF(J.EQ.N) GO TO 310 NEW 

JN-J+1 NEU 

DO 311 I-JN,H NEW 

IFCl-J.GT.IL(I)-l) GO TO 311 NEW 

rC-IB(I)+IL(I)-l+J-I NEW 

• Zfl)'"Z(I)-Z(J)*AR(K) NEW 

AR(K)-AR(K)+BETA*Z(I) NEW 

311 CONTINUE NEW 

310 CONTINUE NEW 

Z(l;N)-0.0 NEW 

IFCNFACT.EQ.3) GO TO 320 NEW 

NFACT-1 NEW 

GO TO 160 NEW 

320 CONTINUE NEW 

NFACT-1 NEW 

RETURN NEW 

999 IERR-1 

RETURN 

END QSOLV 


SUBROUTINE ZEROLV(A,L) 
DIMENSION A(L) 

DATA LPAGE/ 65535/ 
IF(L.LE. LPAGE) GO TO 10 
N-L/ LPAGE 

LEFT“L-(L/LPAGE)*LPAGE 
DO 20 I-1,N 
LFIRST-LPAGE*( I-l)+l 
A(LFIRST; LPAGE) -0.0 

20 CONTINUE 

LFIRST«LPAGE*N+1 

A(LFIRST;LEFT)-0.0 

RETURN 

10 A(l;L)-0.0 

RETURN 
END 


55 



REFERENCES 


1. 0. C. Zienkiewicz, "The Finite Element Method in Engineering Science," 
McGraw-Hill Book Company, New York, 1971. 

2. N. Levy and P. V. Marcal, "Three-Dimensional Elastic-Plastic Stress and 
Strain Analysis for Fracture Mechanics," Division of Engineering, Brown 
University, HSST Technical Report No. 12, Dec. 1970. 

3. C. S. Desai and J. F. Abel, "Introduction to the Finite Element Method, 
von Nostrand Reinhold Company, New York, 1972. 

4. G. G. Pope, "A Discrete Element Method for Analysis of Plane El asto- 
Plastic Strain Problems," R, A. F. Farnborough T. R. 65028 (1965). 

5. T. L. Swedlow, M. L. Williams and W. M. Yang, "El asto-Pl astic Stresses 
in Cracked Plates," Calcit, Report SM, 65-19, California Institute of 
Technology (1965) . 

6. P. V. Marcal and I. P. King, "Elastic-Plastic Analysis of Two 

Dimensional Stress Systems by the Finite Method," Int. J. Mech. Sc., 9, 
143-155 (1967). 

7. S. F. Reyes and D. U. Deere, "El asto-Pl astic Analysis of Underground 

Openings by the Finite Element Method," Proc. 1st Ins. Congr. Rock 

Mechanics, 11, 477-86, Lisbon (1966). 

8. E. P. Popov, M. Khojasteh-Bakht and S. Yaghmai, "Bending of Circular 

Plates of Hardening Material," Intern. J. Sol. Struc., 3, 975-988 
(1967). 


56 



